{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. 以Lena为原始图像，通过OpenCV实现平均滤波，高斯滤波及中值滤波，比较滤波结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2 as cv\n",
    "filename = \"lena.jpg\"\n",
    "img = cv.imread(filename)\n",
    "average = cv.blur(img, (3,3))\n",
    "gaussian = cv.GaussianBlur(img, (3,3), 0)\n",
    "median = cv.medianBlur(img, 3)\n",
    "\n",
    "cv.imshow(\"origin\", img)\n",
    "cv.imshow(\"average\", average)\n",
    "cv.imshow(\"gaussian\", gaussian)\n",
    "cv.imshow(\"median\", median)\n",
    "\n",
    "cv.imwrite(\"lena_average.jpg\", average)\n",
    "cv.imwrite(\"lena_gaussian.jpg\", gaussian)\n",
    "cv.imwrite(\"lena_median.jpg\", median)\n",
    "\n",
    "key = cv.waitKey()\n",
    "while (key != 32):\n",
    "    key = cv.waitKey()\n",
    "cv.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "以下滤波全部采用3x3的模板：\n",
    "\n",
    "<table width=\"100%\">\n",
    "    <tr>\n",
    "        <th width=\"50%\" align=\"center\">原始图像</th>\n",
    "        <th width=\"50%\" align=\"center\">平均滤波</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"lena.jpg\" /></td>\n",
    "        <td><img src=\"lena_average.jpg\" /></td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <th width=\"50%\" align=\"center\">高斯滤波</th>\n",
    "        <th width=\"50%\" align=\"center\">中值滤波</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"lena_gaussian.jpg\" /></td>\n",
    "        <td><img src=\"lena_median.jpg\" /></td>\n",
    "    </tr>\n",
    "</table>\n",
    "\n",
    "从滤波的结果看，平均滤波后的图像平滑但最模糊，高斯滤波得到了比较好的结果，中值滤波对椒盐噪声的抑制比较好，但应用于lena图像，相比高斯滤波效果并不明显，还丢失了一些细节，例如发梢的头发、瞳孔的光斑，帽子上的毛线等；"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2. 以Lena为原始图像，通过OpenCV使用Sobel及Canny算子检测，比较边缘检测结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "\n",
    "filename = \"lena.jpg\"\n",
    "origin = cv.imread(filename)\n",
    "gray = cv.cvtColor(origin, cv.COLOR_BGR2GRAY)\n",
    "cv.imwrite(\"gray.jpg\", gray);\n",
    "\n",
    "img = origin\n",
    "sobelX = cv.Sobel(img, cv.CV_8U, 1, 0, ksize=3, scale=0.4, delta=128);\n",
    "sobelY = cv.Sobel(img, cv.CV_8U, 0, 1, ksize=3, scale=0.4, delta=128);\n",
    "sobelXY = cv.Sobel(img, cv.CV_8U, 1, 1, ksize=3, scale=1.5, delta=128);\n",
    "sobelX2 = cv.Sobel(img, cv.CV_8U, 2, 0, ksize=3, scale=0.4, delta=128);\n",
    "sobelY2 = cv.Sobel(img, cv.CV_8U, 0, 2, ksize=3, scale=0.4, delta=128);\n",
    "sobelXY2 = cv.Sobel(img, cv.CV_8U, 2, 2, ksize=3, scale=1.5, delta=128);\n",
    "\n",
    "cv.imwrite(\"sobelX.jpg\", sobelX)\n",
    "cv.imwrite(\"sobelY.jpg\", sobelY)\n",
    "cv.imwrite(\"sobelXY.jpg\", sobelXY)\n",
    "cv.imwrite(\"sobelX2.jpg\", sobelX2)\n",
    "cv.imwrite(\"sobelY2.jpg\", sobelY2)\n",
    "cv.imwrite(\"sobelXY2.jpg\", sobelXY2)\n",
    "\n",
    "cv.imshow(\"origin\", origin)\n",
    "cv.imshow(\"sobelX\", sobelX)\n",
    "cv.imshow(\"sobelY\", sobelY)\n",
    "cv.imshow(\"sobelXY\", sobelXY)\n",
    "cv.imshow(\"sobelX2\", sobelX2)\n",
    "cv.imshow(\"sobelY2\", sobelY2)\n",
    "cv.imshow(\"sobelXY2\", sobelXY2)\n",
    "\n",
    "img = gray\n",
    "sobelGrayX = cv.Sobel(img, cv.CV_8U, 1, 0, ksize=3, scale=0.4, delta=128);\n",
    "sobelGrayY = cv.Sobel(img, cv.CV_8U, 0, 1, ksize=3, scale=0.4, delta=128);\n",
    "sobelGrayXY = cv.Sobel(img, cv.CV_8U, 1, 1, ksize=3, scale=1.5, delta=128);\n",
    "sobelGrayX2 = cv.Sobel(img, cv.CV_8U, 2, 0, ksize=3, scale=0.4, delta=128);\n",
    "sobelGrayY2 = cv.Sobel(img, cv.CV_8U, 0, 2, ksize=3, scale=0.4, delta=128);\n",
    "sobelGrayXY2 = cv.Sobel(img, cv.CV_8U, 2, 2, ksize=3, scale=1.5, delta=128);\n",
    "\n",
    "cv.imwrite(\"sobelGrayX.jpg\", sobelGrayX)\n",
    "cv.imwrite(\"sobelGrayY.jpg\", sobelGrayY)\n",
    "cv.imwrite(\"sobelGrayXY.jpg\", sobelGrayXY)\n",
    "cv.imwrite(\"sobelGrayX2.jpg\", sobelGrayX2)\n",
    "cv.imwrite(\"sobelGrayY2.jpg\", sobelGrayY2)\n",
    "cv.imwrite(\"sobelGrayXY2.jpg\", sobelGrayXY2)\n",
    "\n",
    "cv.imshow(\"gray\", gray)\n",
    "cv.imshow(\"sobelGrayX\", sobelGrayX)\n",
    "cv.imshow(\"sobelGrayY\", sobelGrayY)\n",
    "cv.imshow(\"sobelGrayXY\", sobelGrayXY)\n",
    "cv.imshow(\"sobelGrayX2\", sobelGrayX2)\n",
    "cv.imshow(\"sobelGrayY2\", sobelGrayY2)\n",
    "cv.imshow(\"sobelGrayXY2\", sobelGrayXY2)\n",
    "\n",
    "key = cv.waitKey()\n",
    "while (key != 32):\n",
    "    key = cv.waitKey()\n",
    "cv.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "\n",
    "filename = \"lena.jpg\"\n",
    "origin = cv.imread(filename)\n",
    "gray = cv.cvtColor(origin, cv.COLOR_BGR2GRAY)\n",
    "cv.imwrite(\"gray.jpg\", gray);\n",
    "\n",
    "sobelX = cv.Sobel(gray, cv.CV_16S, 1, 0, ksize=3)\n",
    "sobelY = cv.Sobel(gray, cv.CV_16S, 0, 1, ksize=3)\n",
    "sobel = np.abs(sobelX) + np.abs(sobelY)\n",
    "\n",
    "sobmin, sobmax, locMin, locMax = cv.minMaxLoc(sobel)\n",
    "sobelImage = cv.convertScaleAbs(sobel, alpha=-255.0/sobmax, beta=255)\n",
    "cv.imwrite(\"sobelImage.jpg\", sobelImage);\n",
    "cv.imshow(\"sobelImage\", sobelImage)\n",
    "\n",
    "key = cv.waitKey()\n",
    "while (key != 32):\n",
    "    key = cv.waitKey()\n",
    "cv.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "\n",
    "filename = \"lena.jpg\"\n",
    "origin = cv.imread(filename)\n",
    "gray = cv.cvtColor(origin, cv.COLOR_BGR2GRAY)\n",
    "cv.imwrite(\"gray.jpg\", gray);\n",
    "\n",
    "img = origin\n",
    "canny_200_300 = cv.Canny(img, 200, 300);\n",
    "canny_100_300 = cv.Canny(img, 100, 300);\n",
    "canny_200_400 = cv.Canny(img, 200, 400);\n",
    "\n",
    "img = gray\n",
    "canny_gray_200_300 = cv.Canny(img, 200, 300);\n",
    "canny_gray_100_300 = cv.Canny(img, 100, 300);\n",
    "canny_gray_200_400 = cv.Canny(img, 200, 400);\n",
    "\n",
    "cv.imwrite(\"canny_200_300.jpg\", canny_200_300);\n",
    "cv.imwrite(\"canny_100_300.jpg\", canny_100_300);\n",
    "cv.imwrite(\"canny_200_400.jpg\", canny_200_400);\n",
    "cv.imwrite(\"canny_gray_200_300.jpg\", canny_gray_200_300);\n",
    "cv.imwrite(\"canny_gray_100_300.jpg\", canny_gray_100_300);\n",
    "cv.imwrite(\"canny_gray_200_400.jpg\", canny_gray_200_400);\n",
    "\n",
    "cv.imshow(\"origin\", origin);\n",
    "cv.imshow(\"gray\", gray);\n",
    "cv.imshow(\"canny_200_300\", canny_200_300);\n",
    "cv.imshow(\"canny_100_300\", canny_100_300);\n",
    "cv.imshow(\"canny_200_400\", canny_200_400);\n",
    "cv.imshow(\"canny_gray_200_300\", canny_gray_200_300);\n",
    "cv.imshow(\"canny_gray_100_300\", canny_gray_100_300);\n",
    "cv.imshow(\"canny_gray_200_400\", canny_gray_200_400);\n",
    "\n",
    "key = cv.waitKey()\n",
    "while (key != 32):\n",
    "    key = cv.waitKey()\n",
    "cv.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table width=\"100%\">\n",
    "    <tr>\n",
    "        <th width=\"50%\" align=\"center\">彩色图像</th>\n",
    "        <th width=\"50%\" align=\"center\">灰度图像</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"lena.jpg\" /></td>\n",
    "        <td><img src=\"gray.jpg\" /></td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <th width=\"50%\" align=\"center\">Sobel边缘检测</th>\n",
    "        <th width=\"50%\" align=\"center\">Canny边缘检测($\\tau_1=200$,$\\tau_2=300$)</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"sobelImage.jpg\" /></td>\n",
    "        <td><img src=\"canny_200_300.jpg\" /></td>\n",
    "    </tr>\n",
    "</table>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "3. 在OpenCV安装目录下找到课程对应演示图片(安装目录\\sources\\samples\\data)，首先计算灰度直方图，进一步使用大津算法进行分割，并比较分析分割结果。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "scrolled": true
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAAD8CAYAAAARze3ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2df5Bc1XXnv2e6BwkJAQLCDwtlQUKwIYmNHa1MgATbYCOIY+HCG3DFhFQlRcrBVY5hqxZ2U1jeFGV7a5OpuOx1CgobGVgjsHGs2GCCBC4igySwjWQwSB4LLAQqyQgLSSNmmJ65+8f0HZ0+c869b368fj2j86nq6tfv3R/n3tfv+8695/2gEAIcx3Ecm66qDXAcx+l0XCgdx3EyuFA6juNkcKF0HMfJ4ELpOI6TwYXScRwnQ2lCSUTLiWgrEfUS0c1l1eM4jlM2VMZ1lERUA7ANwAcB7ATwNICPhxB+PuWVOY7jlExZHuUyAL0hhO0hhLcB3AdgRUl1OY7jlEq9pHIXAHiF/d4J4L1W4mOOOSbMnz+/JFMcx3Hy7Ny58/UQwm9p28oSSlLWtYzxieh6ANcDwPz583HTTTeVZIrjOE6ez3zmM7+ytpU19N4JYCH7fTqA13iCEMLtIYSlIYSlc+fOLckMx3GcyVOWUD4NYAkRnUlERwG4BsCakupyHMcplVKG3iGEBhF9CsAjAGoAvhZCeL6MuhzHccqmrDlKhBAeAvBQWeU7juO0C78zx3EcJ4MLpeM4TgYXSsdxnAwulI7jOBlcKB3HcTK4UDqO42RwoXQcx8ngQuk4jpPBhdJxHCeDC6XjONOKRqOB1atXt7VOF0rHcaYVn/vc51Cr1dpaZ2n3ejuO40wlP/zhD7FmzRqEEPCxj32srXW7UDqO0/GsWrUKzz77bGX1+9DbcZyO5sUXX6xUJAH3KB3H6VBuuukmDA8Pq9sOHDiAefPmtc0W9ygdx+k4Hn74YVMkAeDBBx9sozXuUTqO02E88MADePLJJ5NptmzZ0iZrRnChdBynY/jKV76C3t7ebLqUt1kGLpSO41TO3r178eUvfxn79u2r2hQVF0rHcSqj0WjgzjvvxIsvvli1KUlcKB3HqYSnnnoK3//+99HX11e1KVlcKB3HaTuNRgP3339/1WYUxoXScZy2cuONNyKEULUZ48Kvo3Qcp23s2bNn2okk4ELpOE6b6O/vx+233161GRPCh96O45TOZz7zmapNmBTuUTqOUyp33HFH1SZMmkl5lET0MoADAIYANEIIS4noBACrAZwB4GUAfxZC+M3kzHQcZ7px66234sCBA1WbMSVMhUf5/hDCeSGEpc3fNwNYF0JYAmBd87fjOEcQq1atmjEiCZQz9F4BYFVzeRWAK0uow3GcDqXqh+yWwWSFMgD4dyL6MRFd31x3SghhFwA0v0+eZB2O40wT2imS3/rWt9pSDzB5obwwhPAeAJcDuIGI/rhoRiK6noieIaJnpsMtTI7jpNm8eTM2b97ctvo2btzYtromJZQhhNea33sAfAfAMgC7ieg0AGh+7zHy3h5CWBpCWDp37tzJmOE4TsX09/fjwQcfbOvF5I1Go211TTjqTURzAXSFEA40lz8E4H8BWAPgOgBfaH5/dyoMdRynM/n617/e9gfptpvJXB50CoDvEFEs5/+FEH5ARE8DuJ+I/grADgD/dfJmOo7TiXz+85/Hnj3qoHFGMWGhDCFsB/AuZf1eAJdMxijHcTqXTZs24eGHH+6Ih+wODAxg1qxZpdfjd+Y4jjMuvvnNb3aESALAc88915Z6XCgdxynE0NBQx92O+B//8R9tqccfiuE4TpZOfajFr371q7bU4x6l4zhJVq1alU80w3GhdBzHZCbejjgRXCgdx1HZsmXLtBDJ733ve6XX4ULpOM4YhoaGsGbNmqrNKEQ7AjoezHEcp4X77ruvrfdRT5a333679Drco3Qcp4XpJJLtwoXScZxR7rnnnqpN6Eh86O04Du69914888wzVZvRsbhQJrjxxhsBAETU8unq6mr5jh8AqNVqo3n4clfXiPMe88S0Wn5edswf88aPlidSr9db6uL5o018u7Qt1iFtj22S/WG1h/cXtz8u5+qU/a21kfdzzC/7hG/X9qHcxj+8HGmflUfu+1w/pPLGdGWLmItkGhfKBFyo5AEe1/MDWKbXhE4KjxQyrQ5eliZ+nCggvN6YjqeVYm/ZrdUl1wOtAs/TyPrkep5H2my1Q/YrX+b1S1vltlQ91j60xNs6eWh1Wn1s/Z/KxD3JYvgcZQLrIaTy4I7r5AEp18v0HO2A0vIUKWs85aXaKMsrUncszxI8LiBaP8m6NM9Wa6fVBr49hKCKkGVH7pvnkXllnda+1MqzyioDF8liuFAmsLyDVDrtoONiZB2U8bdc1kSDlyW9Og3Lm8y1LXWgat6k1iY+XaCVoQnw8PBwsl+0frPK4vmtk5NWrnYi1MrICVosJ/4H5O+Ud63ZOZXce++9pZRbBf39/aWW70PvBNpwSJt7ir+1fIA9XM4d7NYwkueRNmnl50RS2qnl1/pCs0Xz1FJCZfWVXE71k2YDT8PXyT6Rc6NWu2QerS+0vpP9aPUxt1fbX1PJypUr8eabb05pmVWzcuVKfOELXyitfBfKBNYBw9elhlAhhKxYAMU8Bs2bjOUNDQ2NCVjwdJaHJYeFqWWt7XLZEkvLC+Q2FGkvz2PZEn9rImnZLofk1tykZaslZFKYrf+C9r+xbJ8MAwMD+Ld/+7cZJ5LASNvKxIUyQ04QUgeO5bHJcuPvlJfEkSIio95WvkgUcM3DSrVL2ybr1rynWCdvpxXkkHVoZWtXC/C8lpCm9pfW79p+tAJaqf+JZkvKNm2fTJYnnngCjz76KA4ePDhlZR5JuFAmsA4W66CS4qCJgEwff1siKYVFHqjy8hVZpxzGaXNv1lAyFUEGoF46w+vSRCA1HLaETAvo5LxWS8Ss4bO0b7wiqZ0AtJOGVqe2/6351vGwf/9+fPazn51wfucwLpQFsA5QSwj5NuuA1dJyct5nzt4i5Vk2FvGEeOBIbkvNq1n1SEHX7NTK0LDaKoNdcigt7U6dKDV7tHyybSlPNGXHeNi9eze+9KUv4dChQ+PK59h41DuB5uHJ9fyPr3kjmsdmlRmRgQYrPb92UfNIeHBAlqEFDlJemuVFaeVqnjJHel0pj1HbJm219oFmI29btDUlkFo7ZH1aOdJ+ue94ubVazZxLjtvGwymnnILbbrsNPT09uOiiizBv3rxx5Z+uPPbYY6WV7UKZoIh3AaQPBmvujx/Q1hBXonm1lr2piHGundo3X5ZR9iJeVrTJ8sy1sqw5S3ky4OVbbbP6LIQw5lKdIt5e3BZCSF7OpNkv11k2TwVXXXUVPvvZz6K7u3vKyuxUHn300dLKdqEsQMqrTP2pUwePldYawufK1DxSvs0SGstb5AER2c6hoaFCbdLaqPVbHAJ3dXW1LKds18rM9QNfJ8u0TkK5bantsj7rJKXZPJ7/To5arYaVK1di6dKlkyqn0ynzWkoXygx8WKQd6NoEPx+KxTQ8fRwya0Ou1PCOp9PuXwZGbmGUQ0tNaOQByr1FOZTkn1qthnq9PkZsI1oe6flZnm/qxJMarmt9HMVes4Gv49/SZiuNVm/qRMTbLcu02mSdVCbKnDlz8Od//ufo6enB+9//fhx77LGTLvNIwoUygealaAeCNhyTZcTl4eHh7K10vEyt3hTDw8Pjatd4D0jNbm6/Jia5si3vyRJz/lsGZCzBssrIed+Afd2mZp9Vp9XH2m9LaKeKj3zkI7j11ltx0UUXTXnZMxUXyoKkDtrcgWT9+S3Pypq8j/VEj04bplnzdHJOVObVxF3aLqPGst5U4Ea2MbfN8lI1G/k6TehSeVLrte3Wia3IiUNus+rh+8K6yH2y1Go1XHXVVTjjjDPGHSw6EvHLgxKkDtC4bImi5Q1owyxZdupATl2Ww4fE2vbYJqsOTdiLCEAqOi/rlCcMra/4lEXKq5L5U31ZdJjL8xcNPmltS02paFMysmx5/3xZfPrTnwYwcs3l3Xffjd7e3tLqms74qSRBSiCts3DqT50K1PCDjHsRlleizdFZ2/kBx4fmUmhSQ1QZ6baQQpBrr1yWtkR7c16zbJO0RbPD8oCtvrBssMrQ7LHazPOkToZlceyxx+KGG27Ahz70obbUN93IepRE9DUAHwawJ4Twe811JwBYDeAMAC8D+LMQwm9oZK/+M4ArABwC8JchhJ+UY3r5pLxFa33KO4p5UkNQzesC0kJliaYsS5ZTxLvKfWtlStstj02m04SL367I69K8Tms/aeKXsk87MVlptPpkHmmD1R+5MtrB5ZdfjssvvxwAsHfvXnz5y1/Gvn372lZ/p1LEo7wLwHKx7mYA60IISwCsa/4GgMsBLGl+rgfw1akxszrkH1UTlfg7zi1Zj9DSDgQ5BJV1R4aGhpJeZlHvTWtTLl+tVksGoLRnWlq2piLgfJ1MZ9VltSPl8VuimBNHXrYlqpoQypOVVZbVrnYKJefEE0/E3//93+Oyyy6rpP5OIiuUIYQnALwhVq8AsKq5vArAlWz9N8IIGwAcT0SnTZWxVaKJntzGBc/ySrT8Vj38d+pZkzyd5l1ZaXNimvoty5PbpaDF3/zi7JT3JKcItGUpnFa/FvHKLLGU28cb+LD6OtWfuXTtpFarYfny5bjkkktanp5/pDHRlp8SQtgFACGEXUR0cnP9AgCvsHQ7m+t2yQKI6HqMeJ2YP3/+BM0oF2t4pHkLKWGxJublu1HkwaiJpOaBWt6kHBbyPCmvSrPdqlPzoLq6utThOE+vTT9wG7SpBllHEU9ce4+PLI/XK9vD+yt1R5KcAtCCR6kTlLRBpquaD3/4w/jwhz8MAFizZg0ef/zxii1qL1MdzNH2rPqugRDC7SGEpSGEpXPnzp1iM6YO7QCLWN6cRD7dOqYZz3MHi3h/VqQ1FZHmv1MemSSKoeW1WvlSnpXVT7wN1v5Iiblm+0T6WwuEFSHe5mi1SevLov+LKvjIRz7SscPxwcHBUsqdqFDupuaQuvm9p7l+J4CFLN3pAF6buHnVIv+0KaGz5q00T46XX+SdNdayNo+neUmp8lPei4y+yjTatZdxaK31neVVWfN6qTqt7alodcr7zHly0tvU6rRslx6lrJv/1vp0vMP9drB8+XJceeWV+YRt5sc//nEp5U50D6wBcF1z+ToA32Xr/4JGOB/Am3GIPh2xDua4LXWAW3nieisKLg8obdgtbZTDPmkH97Sk56V5eNYBrrVrPB6ldeCnbIp9wNNokXu+nvcbL9ey25r2kEPpVH7rRCjFmZfFf2t9IdvUaVx88cXo6elBT09P1aaMsnr16lLKzQolEX0TwFMAziGinUT0VwC+AOCDRPQLAB9s/gaAhwBsB9AL4A4Af1uK1W1COygsccitk9tTwsAPkpRXIcuztqXaZ7XJugunaMTaEltLUFJlS8891RZZh3UikOIfl61bIqMdcj1/v3jqxKHZotlv/V86nT/5kz+p2oRSyQZzQggfNzZdoqQNAG6YrFGdhPQsJJo3kEojvZSUgGh5ZHnSVk3w5EGsHZScVB65XdpoDa95Pq1fZBtku+V+SAmTJkqpoa81XNf2jbaPtRNOTpDH26ZO59JLL8Wll16KdevW4Xvf+16ltmzatAnLli2b0jI7fw9UiOUxcY/PSifXWweM/La2y4M2NRcXsbw3bahYpG5Zb5G2p4SLly8/RT1o2S88aCLL0uY3LZHj7bUELyVgMZ28prZIu3InxE7mkksuwezZsyu1Yc2aNVNe5vTbE21mIhFj7l1ZkfAc8g4UyyYZLJL5c0Iot+XuMc4JRsrjssocD7IsOSyXomjl5banvD7rtyxTE1NrJJDaD9Lu6SiW//AP/4BPfepTldXf19c35WVOv73QRlJDMs0biXlS2+K6lEcnxUrzcOSytE2zReaxDn5uF8+n2czL0YbWRWzll8doNsgTQcr75Ousk5e1Dy2bLY9XC+Zxj187gU20TdOJer2OxYsXo6enBytWrKjEhp07d05peS6UCayDKeUxyd/j9a5klHs8NuaEzxJwKdApkZZ2cxvlK2m1ZWmT1ga+jguSNYyPdqTK17414ZfvOtfK09ZH+7THomn15wQ7V+d04X3ve18llxH96Ec/mtLyXCgTyDmtuAzkvTV54FrDce6pRbSDQ8vH0/Ehm+W9aWgHpuYlaaLS1dWF4eHhFo9Q1qvd925FguV2C+2koF3QHcuy2qitkx6hJmDSm9RsT4mtNVVheaXTXSwvvvjittf505/+dErLc6FMYA2HrKFe/C1FxRqOaweVJXApQZXpUh6h5s0A6blJy/uJ+aQ9fNkajksx4NutNlmiIoUrFYHXvElN4FN9KPPJ/4jcT9rUgtYeIv26yekulADQ09ODj3/845gzZ05b6hsYGJjS8lwoE4zXY9BEEhj7egbrj29dXJ6KlFoHtCW+Vr2WoMg6NUGW6XMRdZkuprG8dOvVGqkgknbCsMTWKkOut8Q9LufuGkrVIdFODtOdZcuWYeXKlVWbMSGO3MeBFCB1gKWERfMYtSg2z5u6A0cKj3VnihyeWx6ZdQDmTgqW2ElvWnp2sl2WoHOsNsq+lWWkxFBb1vojdTJM1at519bdQpqt2pB7pohkpLu7Gz09PXj99dfxxS9+EY1Go2qTCuEeZYKcR5ATmqLl8eXcJSW54Tp/T/VE25Tynvi68UZmpdBZeXICbc3HFvHitD62xNyqP+VNW23STqxFy5hpYgkAJ510Ej7/+c/jyiuvnBaPb3OhLEjKg9S8OJlPkgvgWMNNyx7Nlrg9ZVPqwReWmKWeomMJvIX0eq00ls1am1LziHGb9aSilL3RxtSDmVNinWqb9n+YiQLJqdfruPjii3HzzTfnE1dM50t5hWheSPytDZd4HmuOjG+TwzWtTn4ApSKtlsilAhraUJDnt2wvUmeqL6zgjObl5fpQ6zO5j3KCpom71mZNeFMnIn4lQJHpD6uemc6JJ56Inp4eNBoNbNiwAd/+9rerNmkM7lFmSB0IRbwmKVYpz8mqT9ZliYfcXtS7nKgXpIlJri7ttwzY5NJHtMu3LPu4iGqvtcjVJeuV5WqXR3G0gF6RPtei6DOVer2Oiy66CMcdd1zVpozhyNgDU0BOlOJvmccSRuuASL0TJleXZRffXsRDSQ1H+bI2H6qdPDSv0xI2mSYlyHF96gHIVt/LqQ2rTuuSMEmuH7Qnvsvhv7SRiMYI7JHALbfcgvPPP79qM1rwoXcCebCk0mlD5FSZwNhrDKUg8Ei5JgRSiDUvKjWM0wQxFTmPSO/JitLKPtG8Yu0yKKsNXLS0tPKpR3y7lkeWa3l5uQBPSuhlH2jtl+VY/XakMGvWLFx99dW4+uqrMTg4iHXr1uGRRx6p1KYjaw+ME0uItDQR7SDRRCGFNYzN2WGJlKRWq6mX36TKK1IugOQtfHJZy5NLa3l2qaF/ymPTfstyrRMRMDawlPOorbxaJL7o/2Um093djeXL5Utg2497lAXQRMjyNjWvQkM+lFfzlng5sWxtvlPzkCyPKRXAkfXxOnk/pLwnKU5agEfm1TxbTTSsAI7mScq6rBOCJlCyDTkR1+oF9IebaP1n9a1W9pFKT08PhoeHsWHDBjzwwANtr989ygSpy0H49pQHJdNOlJyg8fVaHi1tqkztgJUnhVwwq0i/WIKRqtcSecsj5OtSD73gWCKpndw0z9USSK1Nqf3oInmYrq4uXHDBBZVcd+lCmcDyMqyDWfPG5DIf9vLy4jI/QCyBkOtynlQqT+ogl+2TZfO0Ka9QK9uam0y1S8Nqe/zw94hr7Uh5b9b+TJ2AZBnaqyWsqQBer1WfA9x66634gz/4g7bW6UKZISWSqSEXxxKRIh6pJcqWp5LyZKwharRP2655TjwP7xc5baDZkmqX9pH9ovWhVQe3zRJ3uS+0fKl+tYRU+69ofagFc7R6nMPMmzcPn/jEJ0YftNEOXCgTpP6gqblJ7WVYOe8g5VXx5dTQukhd2sGrtUfLZ9nH18XHrhHpr16VdVliI9PmBFFro3Vi036nbvuUaa0TR07oLJtTdjlpli1bhk9+8pOl1+NCmaHIH9aKuI53Dk+7lISXp3knlucll7UniOdsKnpZSur95rJsK9Ck2SADQilvMSdGkli29k4bbcitRfOlHVab43eqb7X+k9scnbPPPht/+qd/WmodLpQJLI/BugTIenhtRHt6uSYKlvekzRPKA0w7YLXoa6yHC4MUMNkmKSa1Wi05jxfTpMRc61MrvWaf1od8e6pt2nq5L6xtmtDzftU8RGv/yrK19jppPvCBD6CnpwfveMc7SinfhTJDTpD4Nms5Ck5ElqENmS3h08QyNS+o2cUpOk+qfVt2aAIi08v5uRTRo0uVl0qT+i3XW1MgMq/VB6m6Uic6qyxnfNx4441YuHCh+uCUyeDXUSaQf9iUx6gFA/j2lNBaQzdNOOPvlIBqNmt1S6GyPCperna3kFZ2EYHnv6UXqJWreYSynNQzH3NvprTq4vvAmmbR2pgTRL6s7RsXzfFTq9Vw4403Tnm57lEWQJs3iuuBYl5JquwUKW8ulUZ6sTkPJnpk2hAz1S7Ne0zZmksXQhi9v5mXq11mE9PwNshHwFl2aPvMsjfnfWppxpPHKiN1L7jTXtyjzGDNwcVlzTsZHh5Wn2wt08nytbmq+Ft78IP0NC2KbtcOXC6ecpvVL6lhd6puzdMt6n2HEJK3Zaa83Zw9uTal/h/xW85N8jSpwJ3TGbhQJtBETC5r6fgfPfWKh3q9rg7bItYQToqkto2XZR30mmBYYqIFSqSdMp3sE9kfKZG0gjVaHqvdclmrUxtOa8EWWYd81qSWT26XbbDmLKXAumhWT9anJ6KvEdEeInqOrVtJRK8S0bPNzxVs2y1E1EtEW4nosrIMbyfyj83v9rDSAq2P1sr94fkBZR1k2kEfvb6UVyjrsOYLtbJT3g+gDwstL0sTNk0kNfvjd+oxdJbXbXltwNintfNl69ZVjdxJidvI0+e879Q+ctpHkcmPuwBoj+/oCSGc1/w8BABEdC6AawD8bjPP/yWiqQ0/tRHpyfA5ozjU0wTQ8mh4fstD5fXIMjTxknVEUq+fjReFy23W/KQ1/ZAKXsgyNfG0xE/z4GS+VD/kbJD70sob91NK3GTdFlqgSObT7HaR7AyyQhlCeALAGwXLWwHgvhDCQAjhJQC9AJZNwr7KkSJlDXcB/fmSfJkLpIzwamKoDT818Up5n9oBqUXvuT2WAGnDTU1oLSGVy7xfcnmkfXwbP2GlvFZt31hXBWg2yj6JdWt28jI0obbW5068TjVMJpz2KSLaQiND8/nNdQsAvMLS7GyuGwMRXU9EzxDRM319fZMwozy0A4Rv00QiUvQ6Lm3olxJjvp3n0erWPBTtgNdEVzuYrfq0ddZcpZUnNQS1bEidlFI2ypNWXMefLMTTcPtzos5/54JSqb5wb7KzmKhQfhXAYgDnAdgF4B+b67U9q95EG0K4PYSwNISwdO7cuRM0o1y4F8P/vJpnw5EBHOlpynI1r0dbL0VDy1ev11U7rQBFXCfrjL/j9tRcJG+jJvDavKjmuVlt5ttk3bJe/pFTI7n6rKCLVq9Wpjadop30cvs/1UdONUxIKEMIu0MIQyGEYQB34PDweieAhSzp6QBem5yJ1cH/oDGAo4mGJpqWuFgHvEyj2cAPttRdKHFbfKeNNffI1/E6eR7rpGAdwKmD2hq+FgkIyReJWeJliUtOeK16rfYUHWkUKUvLl1vvtJcJCSURncZ+fhRAjIivAXANEc0iojMBLAGwaXImVov0FjSx4b+1189qHodVj+XFaOmKHHiWnREZ2JHeT0zDhTkn5Jp9Obut+UHpdUVyT/ux2mulKSKQVv9Y6a35Vu1by+N0DtnrKInomwDeB+AkItoJ4LMA3kdE52FkWP0ygL8BgBDC80R0P4CfA2gAuCGEMFSO6eVjHfjyQLbEyDrgtd/cu9I8I+55WZ6TNgyNdcl8fL1sn2xHKuChCZwWcZft5PmtIJQVzIqf1O2K1rC66DSA3HdynVVnqk+1/cHTFbHLqYasUIYQtCdj3plIfxuA2yZjVKeROpD5OikofJsUSevPr4mzdhBJrMuBcvkscSBKP1pMipGsT6aRNmiirPWt7Afeh0Xez62t43dO8W2xvJwwatMnvOyUYFoiaeEi2Rn4TaTjIHUApdKlhtwa/ODJPTVHbpPfMaouhY8HFGJZ2rWFOYG33kGTs016dql8Wv+k+iAVaOPtS9lXVLy0+WqZ1uo/66To4th5+C2MGbSoZYSvjx5dKkggy5Mf60DKDZH5rZCWNynL0WzjNliXGGkertVXViScY83jpfpLm8aQ7bDKKhJZt9qU6hcu/NpJQNpjCaPWV1P9yDBn/LhHmcDyrPi2omVojNcLSwmOVY4m3NKbtMqz2isj0FqemJ5H4HO2SxHm+VN2SE/Zsi/V3iLea84T5ciRgByOyzZZV0O4d9kZuEeZIOXdaZ4TT6d5Pjw/Lzd6KpoYasESy7tJeU+8PK3+mCclyNpwXW5LeU7SBv5sS6sdKW8u5RFKD1TuR1mGFfyxvE2Zjn9b3mK87MrqF55f9qtTLe5RZtAOGk5qWCS9KmuYyb2gXL2WJ5Ty7rTlIh6ydfDy7VpZPG/Ku9PmX3ORfeuEJNvI4QEWWWeqb631mvil8qSmb3KRfaczcI+yIPFPGyOmlrcpD/i4PeXJyedXageVdaBpB5XmQVlt0kQgJ5pWuSnvTUaU5fZYdqqdlt1adF7rDxl4svaVtizrtfpd83Kt9lr7ucg0hdNeXCgTaBP3XCS16wslOc/CGvLJMuSBkxruy7qkjZY3Ke2RbbC8V1lPToxSeVKCl8pjiawWwJHl5tpkzR/Kdmn9ZLWN7ytLXHPRe6d9uFBmyB24GtbkvCxHDruk4FmR0yL1y+3y4NPakvOGUnOTGpP13LQTgJYm1U/aNZMpgZVtSl1xYJ1sZDk5gdZ++9xkZ+FzlOOAz3HVarXkfBd/KpAmiKmDIOdVxt9WeakDU+ZPwcDyceIAAB/hSURBVO9v1/LJIa1mj1aP5qlbEeGi/abV19XVhXq9rtqntSfm0epMnUhS+1frC+0kKW2x7HWqwT3KBPKgsA4y6YHISGjKm7CExcpj3aYo67Y8FUsYpHjJYWHqwNYEg+dNvc9cenXaJxVRzwVKtLYUaVPcbvVlTtSKBses/5BWT6dx8OBB3HPPPQCAT3ziEzjmmGMqtqg8XCgTWF6adqtiJHdwpf74OdFI2Zga7luCkRPy8QR1injNKcGT9efW58QLGBt400RUa5OG1TfWlQya3bl0nSyKnNWrV+PQoUMYGhoC0Ujw6Tvf+Q4++tGPzlixdKEsAD+Acvdz89+aiPAycsIiy5EHvXZgaV6mJsCa3TzPeAI8RebxtL5JCWNKgGXbNKHhdytp5cv0Vr2WB2qlielSQqrVmaq3Ezh48CDWrVuHAwcOjApknNYYHh5GX18fVq1ahWOPPRbXXntt1eZOOT5HmSEVlLEOXEnqThZLQDSRtPJo3qRWh0aRg1I7mON67TmRuQM8JXZyu9Yvmjhb7eXbUnOpWrCIp4v5tbSWHZogW8KZ+l9VzcDAAB555BEcPHhwzAmD//e6urpw8OBBDA4OVmzx1ONCWZCcAFp5cgIU4V6IVadVlhURzgUMZDtkvakDXdZleUtW9F6KnbzG0ZpvtQJoVl9Z4m3ZC+jBJm6vfH1H/OTegZ4apsc2d5o32Wg0sHbtWjQaDdTrIwNQ/t+S3wCwZcuWKk0uBRfKBJaXxs+icXv8E8Xt1iUpmvenBRv4b6scjrXeChxEG2U75Ct2U6LBv7WPNiQvImipbdEmK1/qhGMJN28Db6/VHrn/Zd9ZgRmtHzUB1dpUBYcOHcL3v/999Pf3j/5X4pRG7AO+LqbZtGkT7r777kpsLgsXygLEP4F190yRgyG1TdYl02m3SUbvRRu6yrJ4XZqXZtmgYQl8EaRgyX6RL/cqapNmj1aO5tGNtw2pYJ1cb53kiniM1quN28Vbb72FJ598ctQG62Qav/n/i4iwf/9+HDx4sO12l4ULZQLtj5ybn0oNr6zIcGoeVNaheRspOy1PTVvmbbM8US2/NgfLva+4PvUgYKsemT4nbLn9owlckX2XOxHJ3zkhlOVqrxBJtadM3nrrLTz11FMYGBho8bKjPXEdP/HKQGOtVsPq1auxdevWtttfBh71TpDzDoCxZ3x5kGliKNfLddKbtDwhq14+LLbaJO9Xl+XkvLHUkJvnket5XqtdVrlyWCztSfWHXGett05wHDn01wJtsqzUszm1CLncJ+2it7cXv/zlL1vmW4eGhlCv1zE8PDz6ke2Ol2J1dXWh0WhgeHgYb7/9NtauXYu9e/figgsuaFsbysA9ygLwP0TuAQwpz0VbH/90EvlHlPB8KQ+Uf6eG3daDPiwvSROWFFpezYZUO7T6NI9TCqwWZCnieWpimwq4acupYJa0V5Yx3mmByfLSSy9hx44dYx5Czb/l3La0k6eL/+Gf/vSnbbG/TFwoM/A/jfRaeLAj/o5YHk38LcXB8misg8vyQuSQSLNBHpzyj67dpZPz+uTBZHmTMuChCYX8cFt4Ghlckf1v9aH2WwqaTCv3mybYuXqs9mkntNz881QyMDCAJ554Ajt27BitMwZo4j7jwsf/Y/V6HbVareV98jF9TFer1bBq1appPWfpQjkJ+IFp3QFizT3xA11uk+vkAW2ln+gQTvNiLLS5S62clNfHKfrOnZRo5LxFbR5XE3nL1tTUiiX0VjtSdmpp2iGSzz//vPoqE9kvsh80T1meAOP6Q4cO4YEHHsDu3btLbU9ZuFAmsIajwNg/uXXhtVzHsbyMWLeGfMuf5vXIsosO73Le0ngFQBMVrX755G9Zr/a7iOhbNmsCZb3EzToxWSe1omnlNq0PyxZJAHjuuefw9ttvj3p+vG7Ns5Xee8xj3W3G+6+/vx/r168vvU1l4EKZwPIutGGt5pVYARwpGvKAlkELXg4fEqUOaE1U+HBICqisR7ZL224JoWWH1g+8vdp2q1xuY+rkJJdlGdYUgNYXmqhptub2v5zjk7an/htTxeDgIDZv3jxan5x/5MPmuE6m4ddPWt9xaB7bu2fPnml5QboLZUHijo53ZchhLh9CyoBP6kDg35zcwWEdtKm8qcCRtDVVb8pLk23j1wNqYmXZnxLYiOVZayewlHhZdfA2WYIn26S1VxNyqw+saZOpZHBwEL29vaP/Uy0iL+22phjkiSbeeCHn4Xm6jRs3TrvbHF0oC8A9Ghk0AIoNMXlankcTCi29tINvk96kJdRxmGuVb9UvBUjz3qw81jxcqg8jqSG5th9SNmu2pOZ7ZXmpslLCK/NZomn9b8oQzO3bt2NwcNAUPu5FakLJ+97av/LExNvXaDTws5/9bMrbVSZ+HWUCzQuSQytA//NbwyfrzCw9FP7njPnkMJMj54h4XZrXlRoixvK4HRFLgLW2aeu1AIHlYefeMintswTaEiqZR3pQKRuKiqPWXlmOts9yJ5vxEkLAb37zG7z++uuj89zxO8IfmxaviwwhtJxch4eHWy63Gh4eRr1eR6PRGLWXp49l1mq1lnvkN23ahN///d9Hd3f3pNvWDtyjLID2p5fbrHzaQVqkPnnwaGJg2aDVoXlD2kEfkZ5nyj6tnnjAyXK1uqz8VvqUcKRezGUJ0njaxOu38ub+J9byeP9f4+HVV1/F66+/Plqm9AqjmPHtfJ1cz8vQToD85G+lefrpp6ekbe0g61ES0UIA3wBwKoBhALeHEP6ZiE4AsBrAGQBeBvBnIYTf0EhP/DOAKwAcAvCXIYSflGN+ucizeipYoXlXMY32baWxDiL5m//xLC9E81pSHo60yfLgrPZqomsNnTVbLG9N87y0PuM2pcQrdQLT2pTKw/8TRd47rtWpbbPSjZfXXntt9E6Z7u7uMU8+Ag5fScFPjtGrJKJRj5GnjcvRu4x5uSdaq9VGPUr5WpEQArZs2YLNmzfjwgsvxDvf+c4Jt7EdFPEoGwBuCiH8DoDzAdxAROcCuBnAuhDCEgDrmr8B4HIAS5qf6wF8dcqt7gC0p+pELHFMPVMy9cANbWisob0/3EITYesAt0RJExSePneysEQoHmSarbE/tOsvLSHVTiQ58bFEzRJ+66SWOoHJNqX2wUQYGhoa49XLKQtZd5yj1E5S8uTAy5DL8cMvOZLtip9NmzZNuI3tIiuUIYRd0SMMIRwA8AKABQBWAFjVTLYKwJXN5RUAvhFG2ADgeCI6bcotbwOaJxBJPQRWQwu4WHUWPeg0Twqwo9uWTTm7LXtytqcOfG0qoagwazZIIdKQJxJNnFJ2pYb88uCP6/gcX/ydEmvtvzQRsRwaGsLevXsBtN6txb95m7V6rROMlkfOh3Pb5d1i0o7orXYy4wrmENEZAN4NYCOAU0IIu4ARMSWik5vJFgB4hWXb2Vy3a7LGtpurrrqqahM6knjw82Gc0xnEoE3cN/V6HUNDQwAOi1wUrqGhoZagjnx+QNzPIYSWh2IAh0/GUfjjkJs/QCOuT73NM9px55134uqrr+7Yd+4UDuYQ0TEAvg3g70II+1NJlXVj3Cciup6IniGiZ/r6+oqa4TiOwfDw8Oj91NLLk0+h0obM2nSD9Bqldyo9eW0aJuaL6TSPttFo4F//9V879vrKQkJJRN0YEcl7QwgPNlfvjkPq5vee5vqdABay7KcDeE2WGUK4PYSwNISwdO7cuRO133GcJocOHQJw+ALv1OVZfP6QbysyX8q/+UMyYl3xWw7P+VsAuJ0xb19fX8deX5kVShppxZ0AXggh/BPbtAbAdc3l6wB8l63/CxrhfABvxiG64zhTz/DwMAYGBlouFOe3EcbfMkjDPU3pVcZleQuilYcLY6xL5gMw+pQhbgt/e8DTTz/dkV5lEY/yQgDXAvgAET3b/FwB4AsAPkhEvwDwweZvAHgIwHYAvQDuAPC3U2+24zgARucAtcCcdomYvONGXjvJxdIK5FiepOaN8vXWNZlSTNevX99xYpkN5oQQ1kOfdwSAS5T0AcANk7TLcZwEMWAS4dczyrtg4vWOMsDD388dr7Xk3qcM5oQQRgNDMUAjxZjn4Xf4xDxxOW7nwaSuri4MDQ1h27Zt2LZtG6699lrMmTOnXV2axO/McZxpBn/wBNAaMElFl1OXBfFgjwzqyHRaEMcK/vCypAcp7QZaryd+8MEHR+ddq8aF0nGmCUSEo48+2nyPkhzG8jlDGcyRoqkJIR8OW0KolZETzIg2FOfr+/r68NBDD2HHjh0l9Whx/KEYjtPhzJ49G7Nnzx4dQjcajTFPiIp34UhB4i8Ck9c08vzAiFjFOoCx9/vLl4rxpzs1Go2WoT+/yD4O+aMYxgdoxGsuG43G6DWYXJgbjQb27t2Lhx56CI1GA+985zvxh3/4h5U8SMM9SsfpYObNm4djjjlmzHBVC7zIwAn30rQINaC/qkQbbscyrCG49Fo1bzMu82+eV3tAcExbq9WwZcsW3HXXXZUEelwoHadDqdfro8GM1JDXuvZRRpO1tJow8rQR7WJ069ZEWb8U0HiJEH8hGS+bf8uh/8DAAO699962v6jMh96O02Gceuqpo1HmOGyVw2E+jNWeF8mjyvxJQDHt0NDQ6HCXe27yHd7a3GK0Qz7VP/7mdccyY/k8Eh7Lj6IYo+rcrmh3/B1CwIEDB3DnnXfi/PPPx3vf+94S9sBY3KN0nA7h6KOPxmmnnYZZs2ZlvT0tiAKMjULzYW1Mr5XBh+JyGB/LlXfyFBmWS482puXfEVke9361SPmGDRuwcePGKd0HFu5ROk7F1Ot1nHTSSaMPhODeH3/ikXxeZBQe7XrKuJ4/so4Ha3j5/FpJ+aCMmC+uj2XEJ/5owRxeJ39ohrQl5m80Gi1erAw0xTSxTTzws2HDBjz11FP467/+61IfqOEepeNUyMknn4wlS5bghBNOaJmzsy6z0V6hrM0Fxm2aV6r91uYdrfIA/V04WnlxvXzNCV+O14XyNkbkPKl2Nw8RYfXq1RgYGChlHwEulI5TGfPmzcOpp55qvhNbE524XQZl5AMmpBDKvLIebR3/rd3qmNquibG831wOyXnaeNKI5N6fdODAgdHX75aBD70dp828+93vHg1MaENW/moG7dZB4PCrGvitgtwjkwEe4LDIyGEsDxDxoEkcUsdheqw3Du+7u7sxNDQ0Wg4fznNh02zkUwrxBWU8gMRvd4x2xb6IdsT88VbMp556CsuWLStln7lQOk6bqNfrOO2000bfXRO9rCiCcm6QP0g3iosmMny7NscXRS+mkV4gj67LMvl7b+I2OX8o65Tfsn5eR0zH7YnIC9Wl4PNvWeZU40LpOG1g4cKFOOOMM0a9w3hg86eHRxGJnlORy39iAMa6PEc+iEL+5mlkgCgGYaJw8/JjGgAtd9fw19fKu2247bweeckSJ85f8rt5+AlFvi734MGDpQR1fI7ScUqkXq/j7LPPxllnnZUMjmhBloicS9Qu0rYCNNolPXK+UAZm5GVB1nxkXCfT8wCNDMBY7ebbZdn8t3zeZiT+/sEPfjDhfZXCPUrHKYmjjjoKF154Ycu8oxxy54au2na+jXuFXMii58rFJv6W5cj7xnn5EVkHt01ekiTL5998TlQbgvM6+TY+RcHXAWjx0pcvX15gz4wfF0rHKYFLL710dC6SByDkEDfCAxpyeMoFJfVsSFme9Url+OxJAC2eGX+PN3/VLYAxQ3wZMJKEEEbbz++ykfOsMlAUAzryGZbyARryLp7zzz8f5513XqF9MxF86O04JbB27Vo8/PDDLcNOftcJkL7OUbsMRttmpdHu5uH1WtckRuTwXtav2Z67e0cO+2U9Mr02bK/X62r6MkUScKF0nFI5cOBAUsjknKO1TV4nKbfHNFzc5LyhJcjWHGhc5k/10eY0ebu0eUXtBBC38TR8Wd5KqQl5TP+ud71raneagg+9HadE1q9fjxNPPHH02skIH15y+K2C8sEWWvQ7ot06GIfz8tUMWhSaD435dZzcJqD1dsJYD3+VBBd17bZDPtSOZWvTCdxe/nAMAKNp58yZg8suu6wtr4twj9JxSmbv3r1Yv349gLEPhOBDVu5JaV6YFoHmy0Wiy1pEW+aXLyCTnmJR7zKiDcHjt/YADa1OGUUH0DaRBNyjdJy20N/fj1/84hc4++yzx7wYTD6FnF8Azr0peVF49M6k98nLBVov/ObrNA8zpuVeJfdSI0St1zTKIBW/PlRulx4pX475Y9lxHX/82xVXXIGjjz560vtkPLhQOk6b2L59O+r1Os4880wAaBECfveKvJSoq6ur5ZmOPCLMI+IyD9D6/EoO91J5Hjnsjp4dF1Uu7NrdQvKyJADqOm2OMvaDFlXv7u7G4sWL2y6SgAul47SVbdu2YdGiReatgfKyoCgWUjhSYilFLq6X93bH9bGcmIfnlfdocy9SpudlSsHmc6pyvlL2gVbmkiVLsGTJktEheLtxoXScNvODH/wAZ511FhYvXtwiEECrOMQHTkRB4SLIh7FRLGXghpfDPUx5TWTMy6+d1IbmcugtH5TBPc1arYbBwcGWAEwchgMjYjo4ODjmoRwxCBSH/EcddRSWLVuG2bNnT/2OGAculI5TAb29vajX6/jt3/5tAGMfRiGHt3LYC4x9SK/lkXIPLiKH+MDYO4L4XT+axyfrknZEtGG/tEu7bOiP/uiPMGvWrEn39VTgUW/HqYitW7eq0WIZAY7IKLWMhsvrDmVZ2gXm2nYtX0Tmi+u0bTJ6zn/zCLsso7u7G4sWLeoYkQTco3Scyggh4KWXXsKiRYtaos1ymA0cvn4xipYcTvOouHw2pAyKaOXLiLk27OYeKxc2+TzJ2A6eLy5HZNAIAObMmYNFixaV+kqHiZL1KIloIRE9TkQvENHzRPTp5vqVRPQqET3b/FzB8txCRL1EtJWILiuzAY4zndm+fTsGBgZavLeU5yg9Ps2bi7+1NPK6RQBj3qdt1aOVFe3iNkjkHTXa9ZuzZ8/Gueee25EiCRTzKBsAbgoh/ISI5gH4MRE92tzWE0L4PzwxEZ0L4BoAvwvgHQDWEtHZIYTWR384joNGo4HHH38cxx57LC644IIWzxJovR6RX9PIPUg+58eDITx9LFPe6cMDOtyDBFov9Ym28Og1v2uGz2FGe7Unqce00X4iwqmnnoqTTz655J6eHFmPMoSwK4Twk+byAQAvAFiQyLICwH0hhIEQwksAegGU83x2x5kh7N+/v+VCbEC/cyauT3l7mjcq5wXlsyhlXRxrrpPbyOcyZYCHb+N1dHd3Y8GCBR0vksA4gzlEdAaAdwOIL9P9FBFtIaKvEdH85roFAF5h2XYiLayO42AkuBPnHrUAjwysaEIq82qBFxmE0cTSuo1R2sPLkcLLBbera+yThc455xyceOKJZXbplFFYKInoGADfBvB3IYT9AL4KYDGA8wDsAvCPMamSfcyD64joeiJ6hoie6evrG7fhjjPTeOWVV7B27doWodIi1Vx04rJ82k58iyHfXq/XW54QHn/LcuWcpTWHyb/r9XpLnVYd9XodixcvxjnnnKPOZ3YqhSwlom6MiOS9IYQHASCEsDuEMBRCGAZwBw4Pr3cCWMiynw7gNVlmCOH2EMLSEMLSuXPnTqYNjjOj+NWvfqUOh6WHZ/0G7AdoxN8Ra8idG+6nLu+RHiYvc/78+epTkzqdIlFvAnAngBdCCP/E1p/Gkn0UwHPN5TUAriGiWUR0JoAlADZNncmOM7PZunUr9u3bp0awNWHLDa3jcvTypNilBNmKUlvf8ilIse5Zs2ZhwYIFOO6449rZlVNGkaj3hQCuBfAzInq2ue5/APg4EZ2HkWH1ywD+BgBCCM8T0f0Afo6RiPkNHvF2nOKEELBhwwYsXrwYZ5111phtQOt8oHxCD3A4wh3fvR3XAYevs+TXMvJodPzIZ1XG71i/9tqGWD6PsM+ePbtjL/spSlYoQwjroc87PpTIcxuA2yZhl+Mc8fzyl7/EkiVLWi6tiUgPkl+8LYfR8sG3fEjMn/rDL/3hafnthvzSoug5yse38eH/8ccfP63mIi2mfwscZwazdu1aAK1DZx5wiQETLejC0/FAjgzSRMHjZWmBHJleBnZ4HXPmzMEJJ5wwI0QScKF0nI6m0Whg7969LZHuiDUfmbq8SMunpeV5tKeQpwJLc+fO7aj7tKcCF0rH6XC2b98+uiyFC9CvX9QuDLcCNbmIOq+De4jW7ZBctGcKLpSO0+G88cYbo6++lcNgPgTmQ185VOe/4ytfreG1HH5zcZbXaPLts2bNQnd3d9XdVQoulI4zTXj55ZfHDLE1L9FaJwMtsgw+BNcuQ+LXTFqe60zFhdJxpgkvvvgi9u/fr17InbuuMjU0l8NzXm781gQ4Ih/jNhNxoXScaUIIAT/60Y8wNDTUEv3mUWgZ6dYuGufDZaB13lNGx+WtibzOoaGh0ScVzXRcKB1nmrF+/Xr09/erQ+b4zUWRr5PD8PiyLis4pHmpb731Fg4dOjTmaUczGX/CueNMM/r7+/H444/j7LPPxqJFi1q2RTGz3rXN77iJn3q9Pvr0dP7CME5fXx8OHTp0RHiPGu5ROs40JV42xIfeVgBH/k4Ff7Q5z76+viNWJAEXSseZtjQaDbz99tstgR0pePy6RnkXTu5i9Jj+17/+dfsb12G4UDrONOaxxx7D1q1bxwRx+Hyk9jxJHsCJ27u7u1u2DQ8P49VXXz2i5iItXCgdZ5qzfft29Pb2mpcFaXfzAGNvdYzrgJE5yl27drW5JZ2LC6XjzAC2bdum3prIkfOREW2+8o033nBPkuFC6TgzhIcffnjUs0wNs/k8pbxuMoSAbdu2Yd++fVU3p6NwoXScGURvb++Ye721CHiEe5KNRgO9vb0VWt+5+HWUjjPD2Lp1K84555yWy3niNZT8yePxyecDAwPYvHlzVeZOC1woHWeGsX379tG3HQKHn14OHPYg41PLX3nlFezevbtKc6cFPvR2nBnItm3bRt/mCKAluh0/AwMDfvlPQdyjdJwZys9//nOcfvrpLXOSg4ODePLJJ10cx4kLpePMYJ544glcdNFFGBwcxMaNG9Hf31+1SdMSF0rHmcH09/djx44d2LFjh4vkJHChdJwZzrZt26o2YdrjwRzHcZwMLpSO4zgZXCgdx3EyuFA6juNkcKF0HMfJ4ELpOI6TwYXScRwngwul4zhOBuqEN6sR0a8B9AF4vWpbGCfB7UnRafYAnWeT25Om0+z5TyGE39I2dIRQAgARPRNCWFq1HRG3J02n2QN0nk1uT5pOsyeFD70dx3EyuFA6juNk6CShvL1qAwRuT5pOswfoPJvcnjSdZo9Jx8xROo7jdCqd5FE6juN0JJULJREtJ6KtRNRLRDdXZMPLRPQzInqWiJ5prjuBiB4lol80v+eXbMPXiGgPET3H1qk20AhfavbZFiJ6T5vsWUlErzb76VkiuoJtu6Vpz1YiuqwEexYS0eNE9AIRPU9En26ur6SPEvZU0kdENJuINhHR5qY9n2uuP5OINjb7ZzURHdVcP6v5u7e5/YyptCdj011E9BLro/Oa60v/X0+YEEJlHwA1AL8EsAjAUQA2Azi3AjteBnCSWPe/AdzcXL4ZwBdLtuGPAbwHwHM5GwBcAeBhAATgfAAb22TPSgD/TUl7bnPfzQJwZnOf1qbYntMAvKe5PA/Atma9lfRRwp5K+qjZzmOay90ANjbbfT+Aa5rr/wXAJ5vLfwvgX5rL1wBYXcJ/yLLpLgAfU9KX/r+e6Kdqj3IZgN4QwvYQwtsA7gOwomKbIisArGourwJwZZmVhRCeAPBGQRtWAPhGGGEDgOOJ6LQ22GOxAsB9IYSBEMJLAHoxsm+n0p5dIYSfNJcPAHgBwAJU1EcJeyxK7aNmOw82f3Y3PwHABwB8q7le9k/st28BuITiKxvLt8mi9P/1RKlaKBcAeIX93on0n60sAoB/J6IfE9H1zXWnhBB2ASMHBYCTK7DLsqHKfvtUc1j0NTYd0VZ7msPEd2PEQ6m8j4Q9QEV9REQ1InoWwB4Aj2LEa90XQoivXOR1jtrT3P4mgBOn0h7NphBC7KPbmn3UQ0SzpE2KvZVStVBqZ7AqwvAXhhDeA+ByADcQ0R9XYMN4qKrfvgpgMYDzAOwC8I/ttoeIjgHwbQB/F0LYn0raDpsUeyrroxDCUAjhPACnY8Rb/Z1EnW3pH2kTEf0egFsA/GcA/wXACQD+ezttmghVC+VOAAvZ79MBvNZuI0IIrzW/9wD4Dkb+ZLuj29/83tNuuxI2VNJvIYTdzT/+MIA7cHjo2BZ7iKgbI6J0bwjhwebqyvpIs6fqPmrasA/ADzEyz3c8EcWXCPI6R+1pbj8OxadaJmPT8ua0RQghDAD4Oiroo/FStVA+DWBJMzJ3FEYmlde00wAimktE8+IygA8BeK5px3XNZNcB+G477Wpi2bAGwF80o4TnA3gzDj/LRMwXfRQj/RTtuaYZST0TwBIAm6a4bgJwJ4AXQgj/xDZV0keWPVX1ERH9FhEd31w+GsClGJk3fRzAx5rJZP/EfvsYgMdCM6JSsk0vshMbYWTOlPdR2//Xhag6moSRSNc2jMyn/M8K6l+EkWjkZgDPRxswMl+zDsAvmt8nlGzHNzEyVBvEyJn1rywbMDJE+Uqzz34GYGmb7Lm7Wd8WjPypT2Pp/2fTnq0ALi/BnoswMgzbAuDZ5ueKqvooYU8lfQTgnQB+2qz3OQC3sv/3JowEjx4AMKu5fnbzd29z+6IS9pll02PNPnoOwD04HBkv/X890Y/fmeM4jpOh6qG34zhOx+NC6TiOk8GF0nEcJ4MLpeM4TgYXSsdxnAwulI7jOBlcKB3HcTK4UDqO42T4/9AfQtLgYY1yAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZMAAAEWCAYAAACjYXoKAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO3df5wddX3v8df7bH4QEEggCxeS1EBNq8itGFKI4vUXFgLtbeijpTdWS0q5N1WxtT8tWlssSB9q64+mVWpaIok/QIpaYg1NU4RrrYBZIPJTbpaoZElMVhIiggR2z+f+Md+zmd2c3Zxk9+ycOef9fDzO48x85zsz35lM9nO+8/1+ZxQRmJmZjUel6AKYmVn5OZiYmdm4OZiYmdm4OZiYmdm4OZiYmdm4OZiYmdm4OZiYNYGk6yV9oInb/7GkU5u1fbND5WBiLUnSMkl3S3pG0q40/Q5JKrpszSYpJL1kRNr7JX22Nh8RL4qIrQfZzusl9TWrnGZ5DibWciT9EfC3wF8D/w04EXgbcA4wbZR1uiatgAb4nNtwDibWUiQdC1wFvCMibo6IpyNzX0S8JSL2pXzXS7pW0npJzwBvkPSLku6T9CNJ2yS9P7fdr0r63RH7ul/SRcp8LNWA9qb001OeGZI+Iun7adk3JM1Iy/5Z0g9S+tclvXyM4/olSZslPSXpm5J+bpznaaj2IulCSQ9LelrSE5L+WNJRwK3AyemW2I8lnSxpuqSPS9qePh+XND233XdL2pGW/e8R+znUcz4/rX9pWrZH0tsk/Xw6x09J+vvxnAdrIRHhjz8t8wGWAAPAlIPkux7YS1ZbqQBHAK8H/nua/zlgJ3BRyv/rwN259V8BPElW0zkfuAeYCQh4GXBSyvcJ4A5gDtAFvBqYnpb9NnA0MB34OLB5RPk+kKYXAruAs9M2lgPfq22nzrEF8JIRae8HPlsvD7AD+B9pehawME2/HugbsZ2rgLuAE4Bu4JvA1blz/wPg5cCRwGdG7OdQz/n8tP4/pLznAc8B/5L2Pyedl9cVfd35M/6PaybWamYDP4yIgVpC+iX/lKSfSHptLu8tEfFfEVGNiOci4o6IeCDN3w/cALyulhdYIGlBmv9N4AsR8TzwAllQeCmgiHgkInZIqpAFjHdFxBMRMRgR34xUO4qI1ZHVnPaR/bF/RapZjfR/gE9FxN1pG2uAfcDiMc7DvemYn5L0FHDFGHlfAE6TdExE7ImIe8fI+xbgqojYFRH9wF+mcwFZwP10RDwUEc+mZSMdyjmvuTrl/XfgGeCGtP8ngP8EXjlGea0kHEys1TwJzJY0pZYQEa+OiJlpWf6a3ZZfUdLZkm6X1C9pL1k7y+y0jX3ATcBbU5B4M9kvbyLia8Dfk9VCdkpaJemYtO4RwGMjCympS9IHJT0m6UdkNQ1q+xvhxcAfjQgO84CTxzgPCyNiZu0DfHCMvL8KXAh8X9L/lfSqMfKeDHw/N//9XDlOZvg5HXZ+66WNdc5zduamf1Jn/kVjlNdKwsHEWs2dZL/alzaQd+Qjrz8PrAPmRcSxZLdX8r2/1pD9Mj8XeDYi7hzaUMTKiDiT7BbPzwB/AvyQ7LbMT9fZ92+kMr4JOJbslg4j9lezDbgmHxwi4siIuKGBYzyoiNgUEUvJbh39C1nQhAPPD8B2suBW81MpDbLbZXNzy+bV292I+YOdc+sQDibWUiLiKbLbK5+U9GuSXiSpIukM4KiDrH40sDsinpN0Ftkf/Py27wSqwEdItRKA1CB8tqSpZLdhngMGI6IKrAY+mhqvuyS9KjVYH00W9J4ka1/4qzHK9Y/A29I+JOmo1HB9dMMnZhSSpkl6i6RjI+IF4EfAYFq8Ezh+xK23G4D3SeqWNBv4C6DW5fgm4FJJL5N0ZFp2MGOec+scDibWciLiw8AfAu8ma6DdCXwK+FOyBuPRvAO4StLTZH8Ib6qTZy1Zg/Fnc2nHkP3B30N22+dJ4G/Ssj8GHgA2AbuBD5H9v1mb8j4BPEzWqD3a8fSQtZv8fdpHL/BbYxzHofpN4HvpdtvbgLem/X6HLHhsTbfXTgY+APQA96fjujelERG3AiuB21MZazW3fWPsu5Fzbh1AEX45lnUOSZcAKyLiNUWXpdVJehnwIFmvs4GD5bfO5pqJdYx06+YdwKqiy9KqJP1KunU2i6wW9hUHEmuEg4l1BEnnA/1kt8w+X3BxWtnvkJ2nx8jaXt5ebHGsLHyby8zMxs01EzMzG7cpB8/SXmbPnh3z588vuhhmZqVxzz33/DAiusfK03HBZP78+fT09BRdDDOz0pD0/YPl8W0uMzMbNwcTMzMbNwcTMzMbt6YGE0l/IOkhSQ9KukHSEZJOUfYK1i2SviBpWso7Pc33puXzc9t5T0p/NI0XqKUvSWm9ksZ6RLeZmTVR04KJpDnA7wGLIuJ0spcCLSMbVfuxiFhA9pyiy9IqlwF7IuIlwMdSPiSdltZ7OdnLez6ZHrjXRfbI8AuA04A3p7xmZjbJmn2bawowI72b4kiyR1y/Ebg5LV8DXJSml6Z50vJzJSml3xgR+yLiu2QPoDsrfXojYmt6wdGNNPbYcjMzm2BNCybpLWp/AzxOFkT2kr0a9ancs376yF7dSfreltYdSPmPz6ePWGe09ANIWiGpR1JPf3//+A/OzMyGaeZtrllkNYVTyN7gdhTZLamRas9zqfdCnTiM9AMTI1ZFxKKIWNTdPea4G7NCbN72FA8+sbfoYpgdtmbe5noT8N2I6E8v7fkS8GpgZu6VrHPZ/5a3PtKb3dLyY8neHzGUPmKd0dLNSuearz7MX294tOhimB22ZgaTx4HFko5MbR/nkr1E6Hbg11Ke5cAtaXpdmict/1pkT6FcByxLvb1OARYA3yJ7WdGC1DtsGlkj/bomHo9Z0zw/GAxUq0UXw+ywNe1xKhFxt6Sbyd7kNgDcR/Yeia8CN0r6QEq7Lq1yHfAZSb1kNZJlaTsPSbqJLBANAJdHxCCApHcCG8h6iq2OiIeadTxmTRWBH+BtZdbUZ3NFxJXAlSOSt5L1xBqZ9zng4lG2cw1wTZ309cD68ZfUrFjVgKqjiZWYR8CbtYBqBFXHEisxBxOzFlAN8IvqrMwcTMxaQLjNxErOwcSsBYTbTKzkHEzMWoDbTKzsHEzMWkAwyuMbzErCwcSsBVQj3ABvpeZgYtYC3GZiZedgYtYC3JvLys7BxKwFZCPgiy6F2eFzMDFrAW4zsbJzMDFrAW4zsbJzMDFrAW4zsbJzMDFrAX5qsJWdg4lZC6i6ZmIl52Bi1gI8At7KrmnBRNLPStqc+/xI0u9LOk7SRklb0veslF+SVkrqlXS/pIW5bS1P+bdIWp5LP1PSA2mdlen1wGalExG+zWWl1rRgEhGPRsQZEXEGcCbwLPBl4ArgtohYANyW5gEuIHu/+wJgBXAtgKTjyN7WeDbZGxqvrAWglGdFbr0lzToes2Zym4mV3WTd5joXeCwivg8sBdak9DXARWl6KbA2MncBMyWdBJwPbIyI3RGxB9gILEnLjomIOyProL82ty2zUnGbiZXdZAWTZcANafrEiNgBkL5PSOlzgG25dfpS2ljpfXXSDyBphaQeST39/f3jPBSziReBg4mVWtODiaRpwC8D/3ywrHXS4jDSD0yMWBURiyJiUXd390GKYTb5qm4zsZKbjJrJBcC9EbEzze9Mt6hI37tSeh8wL7feXGD7QdLn1kk3Kx3XTKzsJiOYvJn9t7gA1gG1HlnLgVty6ZekXl2Lgb3pNtgG4DxJs1LD+3nAhrTsaUmLUy+uS3LbMisV9+ayspvSzI1LOhL4BeB3cskfBG6SdBnwOHBxSl8PXAj0kvX8uhQgInZLuhrYlPJdFRG70/TbgeuBGcCt6WNWOn5qsJVdU4NJRDwLHD8i7Umy3l0j8wZw+SjbWQ2srpPeA5w+IYU1K5CfGmxl5xHwZi3AI+Ct7BxMzFqA20ys7BxMzFpANaDqRhMrMQcTsxZQjfBtLis1BxOzFuBxJlZ2DiZmBav14nKbiZWZg4lZwWpNJY4lVmYOJmYFc83E2oGDiVnBXDOxduBgYlawqmsm1gYcTMxahEOJlZmDiVnBXDOxduBgYlYwt5lYO3AwMStYvkbiJwdbWTmYmBUsHz/8eC4rKwcTs4LlayNuN7GyamowkTRT0s2SviPpEUmvknScpI2StqTvWSmvJK2U1CvpfkkLc9tZnvJvkbQ8l36mpAfSOivT63vNSiVfG3EssbJqds3kb4F/i4iXAq8AHgGuAG6LiAXAbWke4AJgQfqsAK4FkHQccCVwNnAWcGUtAKU8K3LrLWny8ZhNONdMrB00LZhIOgZ4LXAdQEQ8HxFPAUuBNSnbGuCiNL0UWBuZu4CZkk4Czgc2RsTuiNgDbASWpGXHRMSd6ZW/a3PbMisN10ysHTSzZnIq0A98WtJ9kv5J0lHAiRGxAyB9n5DyzwG25dbvS2ljpffVST+ApBWSeiT19Pf3j//IzCZQvmbit5pYWTUzmEwBFgLXRsQrgWfYf0urnnrtHXEY6QcmRqyKiEURsai7u3vsUptNsvxF695cVlbNDCZ9QF9E3J3mbyYLLjvTLSrS965c/nm59ecC2w+SPrdOulmpVN1mYm2gacEkIn4AbJP0synpXOBhYB1Q65G1HLglTa8DLkm9uhYDe9NtsA3AeZJmpYb384ANadnTkhanXlyX5LZlVhrD2kyqxZXDbDymNHn7vwt8TtI0YCtwKVkAu0nSZcDjwMUp73rgQqAXeDblJSJ2S7oa2JTyXRURu9P024HrgRnAreljVirVqttMrPyaGkwiYjOwqM6ic+vkDeDyUbazGlhdJ70HOH2cxTRrGW4zsbLyCHizgrnNxNqBg4lZwTzOxNqBg4lZwcJPDbY24GBiVrB8zcRtJlZWDiZmBfMIeGsHDiZmBfMIeGsHDiZmBRvWm8vRxErKwcSsYNXcqHe3v1tZOZiYFazqNhNrAw4mZi3Ed7msrBxMzArmEfDWDhxMzArmEfDWDhxMzArmEfDWDhxMzArmEfDWDhxMzArmEfDWDhxMzAo2bAS837RoJdXUYCLpe5IekLRZUk9KO07SRklb0veslC5JKyX1Srpf0sLcdpan/FskLc+ln5m235vWVTOPx6wZ8qPe3ZvLymoyaiZviIgzIqL2xsUrgNsiYgFwW5oHuABYkD4rgGshCz7AlcDZwFnAlbUAlPKsyK23pPmHYzax3E5i7aCI21xLgTVpeg1wUS59bWTuAmZKOgk4H9gYEbsjYg+wEViSlh0TEXemV/6uzW3LrDTC40ysDTQ7mATw75LukbQipZ0YETsA0vcJKX0OsC23bl9KGyu9r076ASStkNQjqae/v3+ch2Q2sfzUYGsHU5q8/XMiYrukE4CNkr4zRt567R1xGOkHJkasAlYBLFq0yP9draV4BLy1g4PWTCRdLOnoNP0+SV/KN46PJSK2p+9dwJfJ2jx2pltUpO9dKXsfMC+3+lxg+0HS59ZJNysVj4C3dtDIba4/j4inJb2GrP1iDalxfCySjsoFoaOA84AHgXVArUfWcuCWNL0OuCT16loM7E23wTYA50malRrezwM2pGVPS1qcenFdktuWWWl4BLy1g0Zucw2m718Ero2IWyS9v4H1TgS+nHrrTgE+HxH/JmkTcJOky4DHgYtT/vXAhUAv8CxwKUBE7JZ0NbAp5bsqInan6bcD1wMzgFvTx6xUwiPgrQ00EkyekPQp4E3AhyRNp4EaTURsBV5RJ/1J4Nw66QFcPsq2VgOr66T3AKcfrCxmrazqmom1gUZuc/062a2mJRHxFHAc8CdNLZVZB3HNxNrBqDWTNFiw5o5c2j6gp7nFMuscrplYOxjrNtc9jN0F99SmlMiswwzrzVVcMczGZdRgEhGnTGZBzDqVR8BbO2hknIkkvVXSn6f5n5J0VvOLZtYZPALe2kEjDfCfBF4F/Eaafxr4RNNKZNZhPALe2kEjXYPPjoiFku4DiIg9kqY1uVxmHWNYbcSxxEqqkZrJC5K6SJe5pG7Ar/AxmyBuM7F20EgwWUn2XK0TJF0DfAP4q6aWyqyDeJyJtYOD3uaKiM9Juods1LqAiyLikaaXzKxDeJyJtYODBhNJb4qI/wC+k0tbHhFrxljNzBrkmom1g0Zuc/2FpGvTU4BPlPQV4H82u2BmncI1E2sHjQST1wGPAZvJ2ks+HxG/1tRSmXWQ8Ah4awONBJNZwNlkAWUf8OL0/hAzmwAeZ2LtoJFgchdwa0QsAX4eOBn4r6aWyqyDeAS8tYNGBi2+KSIeB4iInwC/J+m1zS2WWedwm4m1g1FrJpJemiZnS1qY/wA/bnQHkrok3SfpX9P8KZLulrRF0hdqo+klTU/zvWn5/Nw23pPSH5V0fi59SUrrlXTFoR26WWvwO+CtHYxVM/lDYAXwkTrLAnhjg/t4F/AIcEya/xDwsYi4UdI/AJeRvVP+MmBPRLxE0rKU739JOg1YBryc7Bbbf0j6mbStTwC/APQBmySti4iHGyyXWWtwm4m1gVFrJhGxIn2/oc6noUAiaS7Zu+P/Kc2LLAjdnLKsAS5K00vTPGn5uSn/UuDGiNgXEd8le0f8WenTGxFbI+J54MaU16xUqh5nYm1grNtcCyTdIulBSTdImnMY2/848G72P8vreOCpiBhI831AbbtzgG0AafnelH8ofcQ6o6XXO5YVknok9fT39x/GYZg1j9tMrB2M1ZtrNfCvwK8C9wJ/dygblvRLwK6IuCefXCdrHGTZoaYfmBixKiIWRcSi7u7uMUptNvnCbSbWBsZqMzk6Iv4xTf+1pHsPcdvnAL8s6ULgCLI2k48DMyVNSbWPucD2lL8PmAf0SZoCHAvszqXX5NcZLd2sNDzOxNrBWDWTIyS9MteDa8aI+TFFxHsiYm5EzCdrQP9aRLwFuB2ojaBfDtySpteledLyr0VW518HLEu9vU4BFgDfAjYBC1LvsGlpH+sO4djNWoJHwFs7GKtmsgP4aG7+B7n5Q+nNNdKfAjdK+gBwH3BdSr8O+IykXrIayTKAiHhI0k3Aw8AAcHlEDAJIeiewAegCVkfEQ4dZJrPCuGZi7WDUYBIRb5ionUTEHcAdaXorWU+skXmeAy4eZf1rgGvqpK8H1k9UOc2K4BHw1g4aeZyKmTVR1S3w1gYcTMwK5veZWDsYa5zJOel7+uQVx6zz+B3w1g7GqpmsTN93TkZBzDqVR8BbOxirN9cLkj4NzJG0cuTCiPi95hXLrHN4BLy1g7GCyS8BbyLrAnzPGPnMbBzc/m7tYKyuwT8kGw/ySER8exLLZNZR3GZi7aCR3lxPSvqypF2Sdkr6YnoasJlNgKpHwFsbaCSYfJrsMSUnkz2V9yspzcwmgEfAWztoJJicEBGfjoiB9Lke8KN3zSZIPnw4llhZNRJM+iW9Nb1+t0vSW4Enm10ws05RjaCrkr1Rwb25rKwaCSa/Dfw62YMed5A90fe3m1kos04SAV3KgonHmVhZjdU1GICIeBz45Ukoi1lHilrNZNBtJlZefjaXWcGqAekul9tMrLQcTMwKVo2gIiG5zcTKy8HErGARIEFFcpuJldZBg4mk9+WmG36CsKQjJH1L0rclPSTpL1P6KZLulrRF0hfSK3dJr+X9gqTetHx+blvvSemPSjo/l74kpfVKuqLRspm1koigUhHCbSZWXmM9gv7dkl7F/ve1w6E9QXgf8MaIeAVwBrBE0mLgQ8DHImIBsAe4LOW/DNgTES8BPpbyIek0slf4vhxYAnyy1k0Z+ARwAXAa8OaU16xUqgEiq5k4lFhZjVUzeZTsNbqnSvpPSauA4yX9bCMbjsyP0+zU9Km9O/7mlL4GuChNL03zpOXnSlJKvzEi9kXEd4Festf+ngX0RsTWiHgeuDHlNSuVfJuJayZWVmMFkz3Ae8n+eL+e/e83uULSNxvZeKpBbAZ2ARuBx4CnImIgZekje0QL6XsbQFq+Fzg+nz5indHS65VjhaQeST39/f2NFN1s0gQgKauZOJZYSY0VTJYAXwV+GvgoWU3gmYi4NCJe3cjGI2IwIs4A5qb1X1YvW/rWKMsONb1eOVZFxKKIWNTd7SfBWGuJCCTcm8tKbdRgEhHvjYhzge8BnyUb4Ngt6RuSvnIoO4mIp4A7gMXATEm1wZJzge1pug+YB5CWHwvszqePWGe0dLNSqVazcSbuzWVl1kjX4A0RsSkiVgF9EfEa4NKDrSSpW9LMND2D7EVbjwC3s79RfzlwS5pel+ZJy78W2c+0dcCy1NvrFGAB8C1gE7Ag9Q6bRtZIv66B4zFrKYHbTKz8Gnmcyrtzs7+V0n7YwLZPAtakXlcV4KaI+FdJD5O9dOsDwH3AdSn/dcBnJPWS1UiWpX09JOkm4GFgALg8IgYBJL0T2AB0Aasj4qEGymXWUmq9uYRHwFt5HTSY5B3KGxcj4n7glXXSt5K1n4xMf46s91i9bV0DXFMnfT2wvtEymbWiakTWAF+R20ystDwC3qxoAZWK20ys3BxMzAo2NM6ErP3ErIwcTMwKNtRm4pqJlZiDiVnBajWTiseZWIk5mJgVLBsBn32q1aJLY3Z4HEzMCha13lyS20ystBxMzArmEfDWDhxMzArmEfDWDhxMzApWq40o6xtsVkoOJmYFi6HeXHLNxErLwcSsYOER8NYGHEzMCjZ8BLxZOTmYmBVs/wh4N8BbeTmYmBWsmh9n4mBiJeVgYtYCKh4BbyXnYGJWsKpHwFsbaFowkTRP0u2SHpH0kKR3pfTjJG2UtCV9z0rpkrRSUq+k+yUtzG1recq/RdLyXPqZkh5I66yUpGYdj1mz1EbA+6nBVmbNrJkMAH8UES8DFgOXSzoNuAK4LSIWALeleYALyN7vvgBYAVwLWfABrgTOJntD45W1AJTyrMitt6SJx2PWFEFWM8le2+toYuXUtGASETsi4t40/TTwCDAHWAqsSdnWABel6aXA2sjcBcyUdBJwPrAxInZHxB5gI7AkLTsmIu6M7H/g2ty2zEqj1purUvE74K28JqXNRNJ8svfB3w2cGBE7IAs4wAkp2xxgW261vpQ2VnpfnfR6+18hqUdST39//3gPx2xCeQS8tYOmBxNJLwK+CPx+RPxorKx10uIw0g9MjFgVEYsiYlF3d/fBimw2qWoj4AVuM7HSamowkTSVLJB8LiK+lJJ3pltUpO9dKb0PmJdbfS6w/SDpc+ukm5XK0Ah4yX25rLSa2ZtLwHXAIxHx0dyidUCtR9Zy4JZc+iWpV9diYG+6DbYBOE/SrNTwfh6wIS17WtLitK9LctsyK41abcSv7bUym9LEbZ8D/CbwgKTNKe29wAeBmyRdBjwOXJyWrQcuBHqBZ4FLASJit6SrgU0p31URsTtNvx24HpgB3Jo+ZqXiNhNrB00LJhHxDeq3awCcWyd/AJePsq3VwOo66T3A6eMoplnhgv0j4B1LrKw8At6sYLUR8HLNxErMwcSsYPvfAe/eXFZeDiZmBQtII+D91GArLwcTs4JFhEfAW+k5mJgVrOreXNYGHEzMClYbAQ9uM7HycjAxK9jw95mYlZODiVnBovbUYI+AtxJzMDErWP7ZXG4zsbJyMDErWG0EfMUj4K3EHEzMCjZ8BHzRpTE7PA4mZgWrVrPncrnNxMrMwcSsBVTSCHi3mVhZOZiYFazqEfDWBhxMzArm3lzWDhxMzAqWfwe8Y4mVVTNf27ta0i5JD+bSjpO0UdKW9D0rpUvSSkm9ku6XtDC3zvKUf4uk5bn0MyU9kNZZmV7da1Y61cAj4K30mlkzuR5YMiLtCuC2iFgA3JbmAS4AFqTPCuBayIIPcCVwNnAWcGUtAKU8K3LrjdyXWSkMPTVY+DaXlVbTgklEfB3YPSJ5KbAmTa8BLsqlr43MXcBMSScB5wMbI2J3ROwBNgJL0rJjIuLO9LrftbltmZWK20ysHUx2m8mJEbEDIH2fkNLnANty+fpS2ljpfXXS65K0QlKPpJ7+/v5xH4TZRPI74K0dtEoDfL32jjiM9LoiYlVELIqIRd3d3YdZRLPmqFZzTw12MLGSmuxgsjPdoiJ970rpfcC8XL65wPaDpM+tk25WOhFZrUS4zcTKa7KDyTqg1iNrOXBLLv2S1KtrMbA33QbbAJwnaVZqeD8P2JCWPS1pcerFdUluW2alkt3mcs3Eym1KszYs6Qbg9cBsSX1kvbI+CNwk6TLgceDilH09cCHQCzwLXAoQEbslXQ1sSvmuiohao/7byXqMzQBuTR+z0smPgHfNxMqqacEkIt48yqJz6+QN4PJRtrMaWF0nvQc4fTxlNGsF1QgqFQF+arCVV6s0wJt1rFqbiZ8abGXmYGJWsIhcm0nRhTE7TA4mZgWregS8tQEHE7OCDRsB70YTKykHE7OCDRsBX3RhzA6Tg4lZgSIiG1vicSZWch0XTJ59frDoIpgNqQWPikfAW8l1XDDZN+BgYq2jFjoqEpWKayZWXh0XTF4Y9P9Wax21mojI2kxcM7Gy6sBgUi26CGZDasGjUhHCNRMrr84LJgMOJtY6asGjNgLeNRMrq84LJr7NZS1kfwO8R8BbuXVgMHHNxFqH20ysXXRcMBmM4NnnB4ouhhmQazNJI+AdS6ysOi6YAOzY+1zRRTAD9ncNrrWZgJ8cbOXUkcHkBw4m1iIi3XWVst5cgN9pYqVU+mAiaYmkRyX1SrqikXVcM7FWsf821/6aidtNrIya9qbFySCpC/gE8AtAH7BJ0rqIeHis9X6w9yeTUbyh5y5VIwjSdzCUln2A3PxgBC8MBi8MVHlhsMrzg9VsfrBKtZptp/a3Joi07sh9BNUqQ9uP2jeRW3f47ZTILa9tj/TSpq6K0j19huUfKkN+fpTtS6KrwtDzp/Jlq+17KC/ZviQNPQCxknZe6/UE8MzzAwwMBlO6xNSu/el5+17IzuG0KRWmT6mkdonhf6xHHsvQ9vcNDt1+qu23qzL8PNSjERlGZs8vfmbfwP7jStFk/QM7srIPVBms1h5Pnz1vpSINveJXB2x5tPI0lK2B7TS4vwby1K6P2hOTp1QqdFWgK30/90KVgWowpZLOeW7/tc4KteKI7NxMrVSYNiV9ui6KOfoAAAd7SURBVCpMn1phaqUydA3X/u260r9jJV3XXblru3ad185z7Tq0gyt1MAHOAnojYiuApBuBpcCowaSrIj719a3c1NM39Mcd9v+hr/1RjtofuhHL0t/vOvlrfxQZtl2zRsw8ciozpnUB8K4bNxdcGhspxfKhYKNcsNkfeLQ/Xy4AZr8R9i/bvz0NTWc5ass0bL/D82jU/ENraYxlh+HU7qMaylf2YDIH2Jab7wPOHplJ0gpgBcDsuadw7ktPGPqlJ4b/+lX6hx+WRv5CGfFrud426uSv/eqsjPprO114FTG9q8LUKWJqV4WpXdmvrCnpl3cq3tCFWNt2/sKt/Yqu7aP2K7Z2Me+/Vkde4PvXr+WpBczBalaTyF/IIy/UYb8Uh+0n285g7pdoRfkayP7/nMMCc6qx5GtXDM3DkdO6mNpV4YXB7FdsrcYRuX3OmNbF1IrYN5DVUPKDBIeXff+xDKTq4lHTpyCUlbsaQ+chr95vhgN/SNSvCdVM6aow//gjATjnJbN5Zt8AXRUN/bsP/UhJ+x6q0Tagkcb8RjbV+I+jBvYX+Vqn0rUVDAym72pwxNQKUypisDq8O3++Bl3bW+1H38Bg8PxAlecHB4dqpAODkfv3Y2i6ts/s33T4Nmt3AIbuJJD/sbj/+qvtd6jmnyvbyLsBtWt36BhGq9lz4DJGHGv+LMfwLENlGr7e+H7Zzpk5o6F8ZQ8m9QLugf+VI1YBqwAWLVoUH1/2ymaXy+ywNPof12wyva+BPGVvgO8D5uXm5wLbCyqLmVnHKnsw2QQskHSKpGnAMmBdwWUyM+s4pb7NFREDkt4JbAC6gNUR8VDBxTIz6zilDiYAEbEeWF90OczMOlnZb3OZmVkLcDAxM7NxczAxM7NxczAxM7NxU6c97lrS08CjRZejYLOBHxZdiIL5HGR8HnwOasY6Dy+OiO6xVi59b67D8GhELCq6EEWS1ONz4HMAPg/gc1Az3vPg21xmZjZuDiZmZjZunRhMVhVdgBbgc+BzUOPz4HNQM67z0HEN8GZmNvE6sWZiZmYTzMHEzMzGrWOCiaQlkh6V1CvpiqLLM1kkfU/SA5I2S+pJacdJ2ihpS/qeVXQ5J5qk1ZJ2SXowl1b3uJVZma6N+yUtLK7kE2eUc/B+SU+k62GzpAtzy96TzsGjks4vptQTS9I8SbdLekTSQ5LeldI77VoY7TxM3PWQvY6yvT9kj6d/DDgVmAZ8Gzit6HJN0rF/D5g9Iu3DwBVp+grgQ0WXswnH/VpgIfDgwY4buBC4lezNnYuBu4sufxPPwfuBP66T97T0/2I6cEr6/9JV9DFMwDk4CViYpo8G/l861k67FkY7DxN2PXRKzeQsoDcitkbE88CNwNKCy1SkpcCaNL0GuKjAsjRFRHwd2D0iebTjXgqsjcxdwExJJ01OSZtnlHMwmqXAjRGxLyK+C/SS/b8ptYjYERH3pumngUeAOXTetTDaeRjNIV8PnRJM5gDbcvN9jH0i20kA/y7pHkkrUtqJEbEDsosMOKGw0k2u0Y67066Pd6ZbOKtztzjb/hxImg+8EribDr4WRpwHmKDroVOCieqkdUqf6HMiYiFwAXC5pNcWXaAW1EnXx7XATwNnADuAj6T0tj4Hkl4EfBH4/Yj40VhZ66S183mYsOuhU4JJHzAvNz8X2F5QWSZVRGxP37uAL5NVVXfWqu7pe1dxJZxUox13x1wfEbEzIgYjogr8I/tvXbTtOZA0lewP6Oci4kspueOuhXrnYSKvh04JJpuABZJOkTQNWAasK7hMTSfpKElH16aB84AHyY59ecq2HLilmBJOutGOex1wSerJsxjYW7sF0m5G3P//FbLrAbJzsEzSdEmnAAuAb012+SaaJAHXAY9ExEdzizrqWhjtPEzo9VB0L4NJ7M1wIVkPhseAPyu6PJN0zKeS9cj4NvBQ7biB44HbgC3p+7iiy9qEY7+BrNr+AtmvrMtGO26yKv0n0rXxALCo6PI38Rx8Jh3j/ekPxkm5/H+WzsGjwAVFl3+CzsFryG7P3A9sTp8LO/BaGO08TNj14MepmJnZuHXKbS4zM2siBxMzMxs3BxMzMxs3BxMzMxs3BxMzMxs3BxOzJpM0mJ7I+m1J90p6dUo/WdLNRZfPbCK4a7BZk0n6cUS8KE2fD7w3Il5XcLHMJpRrJmaT6xhgD2QP3Ku9a0TSb0n6kqR/S+/Y+HBK75J0vaQHlb2X5g8KLLvZqKYUXQCzDjBD0mbgCLL3SrxxlHxnkD3NdR/wqKS/I3ua7ZyIOB1A0sxJKK/ZIXPNxKz5fhIRZ0TES4ElwNr0rKSRbouIvRHxHPAw8GJgK3CqpL+TtAQY64m3ZoVxMDGbRBFxJzAb6K6zeF9uehCYEhF7gFcAdwCXA//U7DKaHQ7f5jKbRJJeSvYa6SeBIxvIPxt4PiK+KOkx4PrmltDs8DiYmDVfrc0EsqfSLo+Iwfp3ug4wB/i0pNpdhPc0o4Bm4+WuwWZmNm5uMzEzs3FzMDEzs3FzMDEzs3FzMDEzs3FzMDEzs3FzMDEzs3FzMDEzs3H7/xGh6IQi1FNVAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "<matplotlib.image.AxesImage at 0x2ade6e20448>"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAUoAAAD8CAYAAAARze3ZAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nO2df5BdVZXvv6u7SXD4Kc+AiInh1/OhovwI6IwMoIhAyhGYGp8wVYkVp4R5QgEz+CPEEGOICA4iM+PTGSlR0HFQGRihiiAQYYhTz5CACYIBjDHRAPJDRSQhJJ3e748+u7N69Vprn9vd957bYX2qbt1z9tk/1t7nnu9Z+8c5l1JKCIIgCGx6mjYgCIKg2wmhDIIgKBBCGQRBUCCEMgiCoEAIZRAEQYEQyiAIggJtE0oiOpWIHiOitUQ0t13lBEEQtBtqxzpKIuoF8DiAkwFsBLACwNkppZ+Ne2FBEARtpl0e5bEA1qaU1qWUtgK4EcDpbSorCIKgrfS1Kd8DAPya7W8E8HYr8mte85o0ffr0NpkSBEFQ5oEHHngupTRFO9YuoSQlbFgfn4jOAXAOAEybNg0rV65skylBEARliGiDdaxdXe+NAKay/dcDeJJHSCl9NaU0I6U0Y8oUVcSDIAi6gnYJ5QoAhxLRgUQ0CcBZAG5tU1lBEARtpS1d75RSPxGdD+AHAHoBXJdSeqQdZQVBELSbdo1RIqV0O4Db25V/EARBp4gnc4IgCAqEUAZBEBQIoQyCICgQQhkEQVAghDIIgqBACGUQBEGBEMogCIICIZRBEAQFQiiDIAgKhFAGQTChuPTSS3H00Ud3tMwQyiAIJhRf/OIX0dPTWekKoQyCYEJw6qmnoqenB5s2bcKKFSs6WnbbXooRBEEwXrzlLW/BI4809wKy8CiDIOhqZs2a1ahIAiGUQRB0Kb29vSAifOtb3xpx7OMf/3hHbQmhDIKg6zjxxBMxMDBgHl+yZEkHrYkxyiAIuoxjjjmm+GeDa9as6ZA1g4RQBkHQNRx44IFYv359MZ7nbbaD6HoHQdA4F110Efbee+9aItkE4VEGQdAYl156KW688UasXbu2aVNcQiiDIGiE97///Vi6dCk2b97ctClFQiiDIOg4l156KW677bamzahNCGUQBB2lp6cHKaWmzWiJmMwJgqBjXHDBBRNOJIEQyiAIOsS8efPUp2wmAtH1DoKg7RBR0yaMifAogyBoK2984xubNmHMjMmjJKL1AP4IYDuA/pTSDCLaB8B3AEwHsB7A/04p/X5sZgZBMNHYY4898OKLLzZtxrgwHh7lu1JKR6SUZlT7cwEsTSkdCmBptR8EwSuIww8/fKcRSaA9Xe/TAVxfbV8P4Iw2lBEEQZdy+OGH4+GHH27ajHFlrEKZANxJRA8Q0TlV2H4ppacAoPred4xlBEEwQeikSB5zzDEdKQcY+6z3O1NKTxLRvgDuIqJH6yashPUcAJg2bdoYzQiCoGk++MEPdvRN5KtWrepYWWPyKFNKT1bfzwC4BcCxAJ4mov0BoPp+xkj71ZTSjJTSjClTpozFjCAIGmbevHlYsmRJRxeT9/f3d6ysUXuURLQbgJ6U0h+r7fcCWATgVgAfAnBF9f398TA0CILu5LDDDsOjj9buTE5IxtL13g/ALdVC0j4A304p3UFEKwB8l4j+BsCvAHxg7GYGQdCNTJkyBc8991zTZrSdUXe9U0rrUkpvqz5vTil9tgr/bUrppJTSodX378bP3CAImuYv//Ivsffee4OIGhfJefPmdaSceIQxCIKWuOWWW5o2YYhOdfnjEcYgCGqxYMGCrnscccWKFR0pJzzKIAiKdOtLLTZu3NiRcsKjDILA5S1veUvTJjROCGUQBCaHH354RxeRdyshlEEQqJx99tkT4pntP//zP297GSGUQRCMYMGCBfjBD37QtBm1uP/++9teRkzmBEEwjCOPPLKjz1GPla1bt7a9jPAogyAYxkQSyU4RQhkEwRBve9vbmjahKwmhDIIARxxxBIgIDz30UNOmdCUxRunQ0zN4HyGiYZ+enp4RYXlBbm9v71Aanj5v57R8W+Yt4/H48ji3RytLi8vrxG3m8TXbvbrXqZMsq9R2Vr5WO2t2eu1QstvKRwu36uz9buqkzcfPOKO9fxSwevXqtuY/0QmP0kH7kfOLXv7o+XddEbHi8TJ6enrQ29trigYv0zum7ed8Zd5aPUr1zgJrpeNtKNNatvL20ISkbv2lbRp16midK+t8auJoleP9JtpF9iQDnxBKB+slpNoPuK6o8HBv34rDRaZOWi66fDvve3WU1Ck75+e1kXbc89Kt/OrA65xSUsVKu/nxONY3T2PZmMvUxF8TfM2OdhKeZD1CKB08T8+Kp13cXIxKnp7cLnlB0qaS4FpiVUeovWOaZ8Xj5a6ylocUYC5olv0lr0wrxxN5z+sr5eG1U24XfgMBgIGBAdMGSbsE88gjj2xLvk3Q7tethVA6aBef1bWUP3hrDEvbt45Jj0OWrR3n+7y77pUjw7RxTKsLKbctsef7JU9S84KttuJjlTJ/bxiBx9HOh6yDdx6t/PkxbYzVaj9v+GI82GuvvUBEO9UyoKuvvrqt+YdQOlgXDA/zvBCtq9eqR9KKjXI/l8+7ecBIL0ezxRIuLa4nhDIvGSd7VqU68rSa3RJNxLT4RCM9fqsdWhmmkGlbaS/5uxkvoZw3bx6OOeYYvPDCC+OSXzfx8ssvtzX/mPUuoF1gpR+xJY6WN5axJiIse7S01j5nYGBAncThZZQufssmzasD7LFLWW7+tiZs+HFPgDX7vPOltbt3/up4oXUEzvptpJTMoZ7RMHPmTNx3333YtGnTuOX5SiI8SgfrQpReltdt4nGsfGU6GQ8Y2Y3mXUZrUsDyDrVJHU9wNHulEFpd+Zxe6+LKWXDZrlY9PPs8m3k4X0mgHS/dFK00rS6fstKXxqbr8LGPfWwojyVLloRIjoHwKGtQ8hzkNt8veTsW2vo9fkyzMeNNnEjbNBstr8tCq7s38ZXrkLct70kTwJzWs0lbfqTZpE0aeWIoy5d18vKx6mS1tffbKnH++efjuuuuw0svvdRSusAmPEoH6wdbx5uQxzwPh38D9oXIL3YuOFwY+ORGPu4Jmdb15XZ4C7FL3f2SZ8Tz0trAE5S6HpsltLwOmt119mX9LIHm5cpjsi29cuvypS99CZs3b0ZKCW9/+9ux++67t5zHROS9731v2/IOoXQoeReel2kdsy5umUYTBs0WjlwXaXWxS/XUbNLKlfZ6y6CsMJmvtr5QE0LNnlK3XSOlNGLpjpdGju3m9NZ59uyxbBqNOFr8+Mc/xt/93d9hl112Gbc8u5Vly5a1Le/oetfAE4fRpPPilrrZdYShJHya9yVtHBgYUNPIsixxtOqoiSDvevOJJstemd4SYK19LBEuea5afqW0PMy7SWl1KdWjFRYtWoStW7fijjvu2KkXmG/ZsqVteYdHWaA0OK9dcLI7Jy9w2fWzBFWb/NDiaMesyR7Pq+Pfvb29I+ouJ25k+2ht5nUpZVtqE00Sq/01m/hTMVo7eOfWmnjyzrtsY0t0vSEaue+tXmiFK664AqtWrUJKCX/2Z3+GPfbYY8x5vpIIoXTQfvDaD1r+8C0vK3eNveNWuaULxhPBUhp5oVpwz68k7jK8ZLu17dUpi6GWlyVYWh5y3xI6rwwvrbUSQuZVp87jwX//93/joosuwrHHHjvuee+shFDWxPsBly7WOh5HHc9BLiDn3o4Wl9sj40vB08Y2Le9Na4M6caw0Wh7ymOd587A69tYNl1i/AS+8zriwDOPnQo5/jheLFi3C8uXLMW3atDF37V8JRAs5aF6AdWFaHy1PQB8z5B6qtrbQmqGWF6Zcc8mxltbwsjSbPA9JE2xNILR0sjsv8+NxtPaSabU20dpbG1bQzm1pVty6Kda5CXjDEvw8tkMoMxs2bMD27dtx8cUXY/r06W0rZ6ITQulQEshSGhle17saGBhwy9PEycpbdoX5I4OWSJTqUdov2aQJr2aTVY4nHJrIW3XQjtcR3ZI9JVHVbLBuOO0USc5VV12FX/7ylzjhhBM6Ut5EoyiURHQdET1DRA+zsH2I6C4i+nn1/eoqnIjon4hoLRE9RERHtdP4dqN5P9qFyD0fKTx11htKT4wvGNcmaTQbeb6WSAIY8TSKZpd18VpLc7R1gLLdrJdA5G1NoDSbtDrK8jXb5XmxbOVtpbWzNVGl3XCsSS+t7bUwaWcnuPfee4fG0S+88ELstddeHSu7m6njUX4DwKkibC6ApSmlQwEsrfYB4DQAh1afcwB8ZXzMbA7Pu5AXVx5bkuvyeHxt+Y91gVgiyfetFzXIvLw6WXF5PO7lynrLF1tYIqAdy3Wqs8yIt3GduFa4J9gc79xY9fFuUp4Ae/Z3Uig511xzDS644AKceOKJjZTfTRSFMqV0H4DfieDTAVxfbV8P4AwWfkMa5McA9iai/cfL2CaxLnx+zPJ6tLzqXuDePhdTTUQtUdIuXK+rqx3XyuJIQcuTEtrkhOVlyeOazdqbf7z0Fp6Y1TneSr51vcROe5MaixYtwj333IPjjjsOfX2v3GXXox2j3C+l9BQAVN/7VuEHAPg1i7exChsBEZ1DRCuJaOWzzz47SjPai9c9KnkLPNwac7L+hsHyVnh6b1xRmyiy1lRKYbc8We8mYOVTqpM3SVOnOy/jldrYO59e/la7aPFzvbwuuvc70urdDSxbtgzbtm1DSgnvfOc7mzan44z3ZI52ZtW+YUrpqymlGSmlGVOmTBlnM8YPz9OxLhyeFtjhXWlLcjzh48j/zJF2tCKqnqhbHpBGSaQ8L1SWB/jvmbTaVh63hNiKa+Gdx7p5ZLTHJHke1hiqVu9u4Ec/+lHXdsfnz5/flnxHK5RPU9Wlrr6fqcI3ApjK4r0ewJOjN69ZtB8t/7ZesKptazPBRPbLYKWHUhIRma50cVli6cWX9Zcz6MDwtZ48XJuEkXlb3p1nj1ePOqKjHbMEvk6Zmm2W96il08Zqu00oAeCee+7BqafKqYvm+elPf9qWfEcrlLcC+FC1/SEA32fhs2mQdwD4Q+6iT0S0uz//8WpdTWub55c9MR7mXai8jJI3aXmElr3e0hw+W63ZUWcm17JJ5ufZJL1W+fcPWtvIOHVmxXnZWhotD+/caW1inS+rTSyvuBtYsmTJiCfNmubWW29tS751lgf9O4D/B+CNRLSRiP4GwBUATiainwM4udoHgNsBrAOwFsC1AD7aFqs7RF1Px7sgZBxtv7TkJ4eV0sk8SmjjlaX8NCEERs70e+3l5WPlXaqzbHteNo+jnU+OtoZV1pHnJZd2eTcLIlK74NJmzfZu56STTmrahLZSnMZKKZ1tHBrRMmnwV3DeWI3qJixviR/X4nK8i6m0BKgkmrIcWb42meF5etLL8jwwSwzlTDzP2xIxWQerrTSbPI9Ns0/WtSSyJdHXyrFusFZ71UnXzdx9990AgJNPPnlouynOPPNM3HLLLeOaZ/f69V2AJYzSmyh5RzzcuijzvuUxamWV1hTWsc+6IOteqF595dt7+DFNNLQZZqssS+z538DKsuX/9oxGZFttF8vT9urlxe127rrrLuy6666N2nDnnXeOe54hlAUsAZH72sWvzQpbZbRathQV+fcPpa6ctTRHS8u3S+JbWmg91otf5iW75doCdq9syzZP3DW8GWuvPeR5qLugvpv52Mc+hjlz5jRW/ubNm8c9zxBKB8/Dsrwi7wLTRIZ3U0vepNzXRLvkIcn0JY9S1s9rCykWngDL4966UNkGddrb6rZrE1o8f6sOWj21m6GsozVcYp07q04TicsuuwzXXXcdUko45ZRTGrHhb//2b8c1vxBKhzoeAd+v48W0cgHUiWtd+Fpe8sK1PBxLPKUt8rjs9mp10OpUWpytiZyVh2Wn983bhA8XWG1ohedlXHXtsFZNlMqbSNxxxx2NLCNasWLFuOb3yn0mqQbaGJvngXkekrZv4c2YW16it3jdu+C8enjeTQ7X1k3yOPK4zMOyyUOzS/6FhFd3L0zWryRgMp/SjcizRSt7ogslMLiMqNP1eOSRR8Y1v/AoHTTPJodrcTJ1JmnkvvSaZF4a2kVZ8gg1r82yWx6ruw5T5uGt1dTaUyvbu2FJ20rx+bc3geNta3ZZNxbrvGp5l4YgJiopJZx55pl41ate1ZHyXn755XHNL4TSwRIa75gmptrbdeqU3eqstWV7nfI0kdQmFqz6e2JXxzuSAsftKv27o2WP5aFZdljHNNHVhF6bVbfqWaKOvRONm2++GX//93/ftBmjIoTSwbrAPO+QH8+eDu8S8idGtFlxT4CsmWxui+f58TWV3ngft0+LKy9c60kZabN24XsTTpbH5rWZdq60+sm0WhorzCrX8py9G4ds37qiPlFZvHgxUkq44IILJtTbiCaOpQ3gXVjWcQ+v656PWwJs7UusJz+0Oml5aReoZq/Mv85F7YmoZqNVljUe63lxMo7X9S8JlSbU2rZWJ224pJTHziaWAPCP//iP2HPPPbFy5Urcfffd6O/vb9okl/Aoa1LH25Db2r6Vr+yO8e56nRldK64UBc8my7OTyD8iK9Xfa4OSyGnHWlmMbglfbl9vXNYqg8/uSxu9G4lXN004d0aB5Fx22WVYsmQJzjuv+x/mC6F00LpIPDxvAyO7itLTKXmEmhho4pyPcXK3Wq69tDwzy5uxhF/aL+tgiQRvr9IEl2aj/JZ5l/5ewruZ8PbyypLnUpbriSIfuvDaS7a1l/fOyDXXXIOUEubPn4/TTjutaXNUoutdwLsQWvFCNI9IS1fqjnoCqqUrzUpbXXiZVpZXGrcsjYFye4jsN597+9aax5KHNjAwMDQ+ZtnvnVP5/+ajETYr7itBGC0uu+wyAMBee+2FF154oWFrhhMeZU0sr6wkBEC5q8jzk2+vsWZIvQu7JO4lYbc8LFmO94qtknh79nqel5bOei+oVnbO33oJb91yJdZfUsjfAD8mZ8k5rxRvUuO8887D0Ucf3bQZwwihdGjFS9AmYrS8rDWUGTkr7nlGsixpb6mr571ZyJut1mbD8z6vj9Um1g1HDl945dcZCijlb8Xlbc2/tbTyXGnngZ8rrU4ynbw5yuf4d3Yuv/xyrFy5EiklfOpTn8KJXfA29eh6O3gXHUcbG7Q8EE1AtXIl1iyttLVVSt5anfQc7QkZKewyTHqEVt48fSluKR/vBuSdO+1GVcrDsrk0jizTvhJZvHgxgOaHJF7ZZ6Em/C4v7/rejLJcQ8nT87TS25FxtUmRnN7zCHl66wKvc2FrnpPlXVmTNlZ+mjh4npbnzXniJ+PKOFrZvA516sGPyfNV8iK1Y1rer1RSSliwYAHe9773NVJ+eJQO+UdqTRpYEyHWD7u0TMe7IOqWbaXR7NcESruwLdu8Y5bNVhxvfWHJbi9c5p8nYkptXaetZHxLLGWeVlzLlmCQz3zmMwCAXXbZpePrLsOjdKjzI89hlteRqTMpw7etFzxo6TVv0RMxyzOzyil5aZY9VltpHptVTp2bgZUu3+T4MetGYOWldZu1+lsTM3I236qTdZ6sY69kLrroIrz1rW/taJkhlAXkBV1HNEpPX0hRlfmWPM/SpAIvWxMlzy5P6LxxUm0tp1Zfz6uqK7IyvrbP01g3MJ7GWu7klWPVRzuv2recyLPOUwjlcP7hH/4Bq1evHnrRRicIoXTwPCdPzAYGBszjnveUj2vixC/okn3WxSnL8i5Cq75aem4TH6ao8zKLusc8T1CzvVXR09ZxWnl6+XpC57V16WYW6Nx8882YPXt228sJoSxQ5wcrRdETO43SC38tD9DzorTtkqfqpfXIaz9LwmAJimVDSTjqlGmR4/Eus7eA3RLSOnXm8WS4tu8dC0Zy/fXX4+STT25rGSGUDpYgyRllHqeOSPK85Fik1iW3vCqZr9fVs+y06sTrrHXJpZ1am8l8PaGxhhNkPE8ctfax1opqH3mu5NuWZP088ZP28W2vl2DVN/C58847kVLCa1/72rbkH0JZoI53o21rF1fe1i5svm+VI9NYwqMd9y42KbDaBIT2Le0seVXcLjk+55Ft8vKTcbQ8tP3SjcVKawmhTCvzLrVVndUMgc1HPvIRvO51rxv3RfohlA7aBeAJkLVekQuDFNxMaZZbW5zM40jPicPXfHpeirY2VKbVjnGbeP1KazllupxGa6fSRJMsR5av2WB5i630GEp18m4emh2y/ereTIJBFi1ahCeeeGLclw/FWaiBdYfXPEbp+VieUClvedy64DQ7LCHhx3ha7pGVvEEtzJq59epmxePPj2uTRNIWKYZae5favc5EjlZvy55W01jllN4tGnSOWHDu4AkOPw7444kyXWkcU/uPbu3FD3UESTte6lrLcjWvxlpy45XtCZZVH8tzlPGknXXSWefGSm9tl+J558vzRoPuIYTSQRM+uc3j8uNWl0xezN44pSUimkDLYzJvLS9vcshajqN9W3X1yvTqp9XLy1OWXxpf5Gm0MUGt7lrbynKtv8SwbpJyWMOqV4hm8xS73kR0HRE9Q0QPs7CFRPQEEa2qPjPZsUuIaC0RPUZEzfz7eZvIP2z5KjRg5HiXTFfH+9KEy7voc751XzPmeVnWH4nV8X68+loiqwmJlpc8ptkpy/bKkMfk29p5HK3r28q5k2Vx+726ynqHUDZPnTHKbwDQ/sH8iymlI6rP7QBARG8CcBaAN1dpvkxEE/YdUdK74WNn8n9btHTWMSkO1gUq8/O8N89+70LzLsiSx6bZVXe80hNpb9LJEj/NbqselgjKfXne5U3JK1M7rtleZ3Y9aJ6iUKaU7gPwu5r5nQ7gxpTSyymlXwJYC+DYMdjXNXAhkD90uV0nrpaWl8PT53D+LcNkfGsmN29bXXbLlrwvZ78t0bBEUvO4vHaw2sSySUujCatmmyXQVrk5nNvszc7XEX9tKCDEsnnGMut9PhE9RINd81dXYQcA+DWLs7EKGwERnUNEK4lo5bPPPjsGM9qLtTzDEotWftTcW+FptW+ZrydUli11RE0TBO27Tr7yWJ2xQ227FVvq2KgJEu8leDcrT+jkvneD8Oy1ygmaY7RC+RUABwM4AsBTAL5QhWtnVl0FnFL6akppRkppxpQpU0ZpRnuRa/Isr8LC8xSsvDTB8i487cKWT5TUWZdo5ZvDZL1k/Xk52trDOmOHXp2t9YtyQkRrc+t8yfapM2ygtQ230Ts/PK98nqw6W20UNMOohDKl9HRKaXtKaQDAtdjRvd4IYCqL+noAT47NxObgP1D+PLOG162WeVnLg7R9LY/s/WhxgB1jfnlNYh2hk+m1C9uapNHy1iiN12rHZJ1aFVitLO88amJaqled9KU86rZD0AyjEkoi2p/tngkgz4jfCuAsIppMRAcCOBTA/WMzsVk0b4Hve96KFdcrR9u38qpz4Xn5Avp/dGviwP9n3BOCUntYdrcyPijt9ijVn9tUEjUZr+5EjNamlgC3KsxBZyiuoySifwdwIoDXENFGAJ8GcCIRHYHBbvV6AOcCQErpESL6LoCfAegHcF5KaXt7TG8/lgh4ax+1fZk272uTNZp3ZNliXeDWMU8IvYkjb5JEy9ebkNDs89ZWajbVPSbL8Oxu9UZXmnzjx2SefL2lVVfLrqAZikKZUjpbCf6aE/+zAD47FqO6DX4h8zGxTEk0eNrSD98SPa0sy07Nk9HESyuPb3uP9ln5yvilNZVavUrib6XXwrUw68/M5GvWLDGW6bw8tPrIfOr+HoLmiGe9R0EWEcB/SqeOMFr55zTegmjP65HppfBpEzyaUPCyNLGV76K06qPlVxKdUn5am2j5e+k9+yy0GyfPyxNJb/w5PMjuJR5hLODNnGqTMtaFYK1r9PKzvCMZx+o2S6+25FVaNnnlW7PZgP5GI83jlPlqbSI/3ppRLS+rra24VlprvFnmZ924tDwtYbXaLOg84VE6lMYiM94P2RM+b19eJK0sTcreoVeOJTIca+KIe9NWHXJaPgNfuvA14fAeWbTytv7oq5VzIcO9G4wWR2t/6x2JnmCGSHYH4VE61PUI+T6Prx3TPDTA977qrAfk4ZZd2oVneTjWBIaXr3WhW7Zn4ajrEZaWF5Xq4pXVav09z9HyFvlf5ZbONce6QQedI85AAc8LseJJ8tIaS7DkukivXOsi0kSwznpMKy0PH83xvG95d9ZaUE/0tfwtUeNkj7Onp2fEMifrhqOVpe3XffZbGx6RNng3nqBZwqOsiealaYKjTeBYXdQcNjAwMKxbZuVd8rjksbyvbWtxtTheWi1/a3LFmg2W35pwyG3ZnnVfOJz/HVMrl8f1/mSs1O51xV3LT5af6xo0Twilg/ZjtbrYlsBo3VUrvRQabod2cdVZ9G5dwJag5WOl45ZXpNkry5O2afuybNlepfpptlpenWZzqb08Mbds8NqoTllBc4RQFvBmL0uiZs0y8/Sed6EtTJZ51NnXxFnG50LtiUCrk0qW/dY2z9MrwxNHyzurK7BWnercECxhH805DG+yewihbBH+lwOyy2xN1HBRHa2HYHkdWn7W2KRlmwWffNDyy9ut/k2FvAHUsTPXqxVRrtNemvB5583yIL36l2zT2mqsv5dgfIlblkO+MK11iPyCL3Wv8j7H8+Isr4/b4XldXlkyrHQxWzZpx7X611nzWPq00uX22tL7SLtKbSmPl0TNK680jt2NfOITn8AhhxyCQw45BJ/4xCeaNqethEfpoHkz3gXPKY3JeeVZi6xLbwKS5VpleherlZdVnnahW+Va3XZLnLS2KtWJYz1pZG3XOTdeO5TOhRWvTvndxIwZM7Bp06ah9k0pYcmSJQCAz3/+8w1b1x7Co3SwLlZrOYrlHcj/7O7p6RkW5l388rg1C87TWGv1LAHW6mytG9XaxZt88G4smjBq4ul1my3htbxbS6jqtFee5JJlaGnz+dW8f629ZXt656gJPvnJT+JP//RP8ba3vQ1btmwZ+g1nW1966SV873vfw+GHH960qW0hhLIGcvlP3rbiap6Nlaa0kN3yzPgxzYOpixQAr06W0EhbShe3J3byeJ38rRtMtj3TyqqNZA8AAB1mSURBVFiqJs6l/9n2bnRamFWHbhHHzLx583Dvvfdiy5Yt6m+OD1Ft3rwZ8+fPb9ji8Se63jXxJkgsgSqJauli0MqpcwFpFyMfZ5Xk9YVaemmHJ2JWeGlsE9D/l1vbtv5x0roJlcTVC5PHsn2W4JYWn9exuzRL3mkuvfRSLFu2DNu3bx9asK8t3M8PEKSU8LOf/axBi9tDeJQFSmONdcfdtDuwLMPzJqWHY63ntNZqSlusOljjlJbHY3U/NWHSBNVLo+Uvu7+lcr16a+FWWOkYz48Pq3hj2bLN63junWTu3Lm4++67sW3btiE7c3c717G3t3fYdk9PD1avXr3TdcFDKGtg/VDrTn5Y3lHpovbgExWWp+kJlWWTV3ZJgEuUhMl6IqZOGZptVtt7N7ISWr48fckz9wSdYz362SkuueQSrFy5csgGzWHQhD2HvfjiizvVTHgIpYO3REfulyY/tDB+EdTpdsofpGWTVqaVpzfZYI2Fyry0i1or33sRME+jtZcVx8qjlDcP1+ptCan3ViZvUqtkjyZEXn3aySWXXIIHHngAL7/88rAJm2yPHMaRN+vsXd52222YPXt2x+1vByGUDvJHIReXWwKkTQSUls9Yd2ig3nsdtfRafbQfvBan1A2s0+UuxdHs9dJa6yGtc1Gyo8550s6nN/ZprQ7w6irTluK1kzlz5mD58uXYtm3biN9dX1/fkI35WpBrifnQyNatW7Fs2TL8xV/8RcfsbxcxmVMT62KQ+1Z31PvRW4Ll2ZJSQm9v79AbxjXqeoSaDV5cebzOhWzl49lgpfXysmarrfpZZVltZtlubddpM02cvTLbxYc//GE88cQTQ5M1+XfGJ3D4RE7e7u3txfbt24fCMvk3+vDDD6vlTSTCo3TQvAGrG8Lj521LPD1hlR6PduHzbf5OR2mfzM+qk8y3lckfy8uy8pcTF9bkU0ncpZ3WudDqZ+1z70ge886BPKdWOVre3hCNZkM7mDdvHmbOnImnnnpqhMeYt/lETl9f35CdPD7/zcoJn8MPP3xCj1mGUI4TlmfRytiTJwqeGIx137sQ5THL+9HstvLgeK8049ueaGg2yuOl/LyJKU/0clpLQOu2gxWnEyL5+OOPqzda66YLYNjYJf/mNvPvzZs347bbbsP555/f1vq0ixDKAnW7g3zbuzC1vFvxHrQX/GppvLxH63nVqb9WtrbN41rvX6zj0dURWK2uEu9/y60bmGQ0AqmJimZ7u3j00UeHLf8BRnr7smfD7ZPto+WR2bp1K1asWNH2OrWDEEoHzcvQ7rjWj8nqjnkiqh3TurWeF9NKt9bzrLRjJeHilOrriaN38yi1V2lbpuPdRO3cemVaQw6WDV5bl4Y8xpP58+fjgx/8IIAd3qHVbeZhvOvN68u73jI/7n0+99xzOPts7R+wu5sQygL8R529H+04/9tWy+uzvB/tQihNPnheTisejLZdSu+VKS92bz2gJuI83LuhyPTakIAluvy4JdJaPjKOvHFY9a0r+lp7toP58+dj/fr1wya56tzU+HEphrk9rN4Kj/Pggw9OuMccY9bboY6oZDTvU/uxaxd2nfjaccvD8R4HlGsZtTqUyrcEhceTwimP15npl3/GJfO3xMSyrc4NxDvmlacdt2ywfiut2jxafvWrX6G/v18tM/9++Kx3Dpez3VY9cx3zjHh+tDGH9ff349FHHx33erWTEEqHuqJodXFzHO2Y5W3wdNrsr9YllPG5zVq5lk08jbTDahPPNq0NtLpbkwB1xsWAkbPVmp2Wzdq2dT4sm62bCW8Tz45SuDVO3goLFy7EH/7wB/z+978fWubDez7bt28fEkYeLvezgGbR44KYlwhx2/kf6/F6rlq1CvPnz8fixYvHXLdOEEJZA0006tzptYug7uRQ6ULn4ZonI8O8CRbtmPcUjZa/tDVfdCUh8dJb8b22l91Jr8xWZ6u1drbCrbhe3qV0Y+E3v/kNtmzZov41BvciczifZOPDCVxg+W9Ziqb0IvM3L3v16tXjUrdOULxVEdFUIrqHiNYQ0SNEdGEVvg8R3UVEP6++X12FExH9ExGtJaKHiOiodleinZS61NZEi+VRSk9B2+flyLLkMTnLqF14pQkcXifrOM9j4cKFqtdjlcFtkGXJ+mjCbdlktb8WX7NHa+fSttzn+XjnQFLnXGi/gVb56Ec/inPPPRf9/f3DJlbktzaBI3/PWniud/6W7+GUZeVPT08P1qxZg4MPPnhCTO7U8en7AVycUjoMwDsAnEdEbwIwF8DSlNKhAJZW+wBwGoBDq885AL4y7lZ3AdrFbh3jYpLxPFN5sWnxtXT8/6u1ONbFqJVpiejChQtHxNEEvs63VnauhzUpUKdcK43WrhaeWJVuaJ7NcnZds1FLNxoWLFigvkPTuyFpi+61m1GOw0Uwh0kBteqZPz/5yU9GXcdOURTKlNJTKaUHq+0/AlgD4AAApwO4vop2PYAzqu3TAdyQBvkxgL2JaP9xt7wDaCdcXjA5rJUfNNHImXF+zBNJ6zj/gZZePqGVaaFNOi1atKhou1YXz3Yepgl2qY1LIg+MXIPqeYCtiKqWNofxMT7enbXyHQ+BBAZF8vnnnx9WT+kR5mNauZpA8nglkZV5aDf6fKy/v39Mde0ELY1REtF0AEcCWA5gv5TSU8CgmBLRvlW0AwD8miXbWIU9NVZjO83y5cubNqErWbBgQdMmBAZ50iaPN/L3AWThtt4PkEW9r68PAwMDIyZ98oQPsGNSJ4fl7zypw/Pg4rh9+/ahMUwedthhh+H9738/rrzyyk40U8vUnk4jot0B/AeAi1JKL3hRlbAR7hMRnUNEK4lo5bPPPlvXjCAIDD796U9j06ZNAOxhCs0btMbhpUepPbaY0bxVnlf+lkMPOX5/fz/uuOOOrl1fWUsoiWgXDIrkv6WUbq6Cn85d6ur7mSp8I4CpLPnrATwp80wpfTWlNCOlNGPKlCmjtT8IgootW7YA2CE+WbDk69K8cceczhpPzPnnb/m+Srmt5ctt5Plt3ry5a9dX1pn1JgBfA7AmpXQ1O3QrgA9V2x8C8H0WPpsGeQeAP+QuehAE48+nP/1pzJs3b9jMshyXlB6lPKZNyPB02iy4FYdox7srtdlyzTPNZa9evborvco6HuU7AcwC8G4iWlV9ZgK4AsDJRPRzACdX+wBwO4B1ANYCuBbAR8ff7CAIgEGR5H8OJ8ceLS9Sm9Thy3s0b9PKK6eXnqw3MSSFk5ezYsWKrhPL4mROSulH0McdAeAkJX4CcN4Y7QqCwGHx4sXDllHxBeN5siSP/eVJHPnETZ5okS/ezfEGBgaGpcv7eVtbSpTjexM3ubyenuGL0LN4rlu3DuvWrUN/fz+uuOIKdAPxUowgmGDkx/6kxyaX5wD+o6nWRy6H07rWcuLHykeObcoX/PLuv7R9yZIlmDt3LrqBeIQxCCYIixYtwuTJk4e8NO6N8YcN5HIg7iUCGBFHhskHF+T7QvmjqdwOTTzztkSGZVHlHufmzZuxdOlSnHvuufjXf/3XcWzJ1gmhDIIu54orrsDkyZOx2267DQmefKggd4vl0pzcFc5ilsczs8hxslBxgcx5cyHN3W6eV06X/yeHd7cHBgbQ19eHlNIwIeSCn8vh6zD7+/vx/PPPY+nSpTjwwANx2GGH4aijjmrkRRohlEHQxVx11VXYbbfdho0PZkGSIig9R+4h8nFBPoYoBTMLIhfGnIZ/88XoPIzbx1+ywfPPIpnHT3OcnAcX0mxjT8/gs+G/+MUvAKDjYhljlEHQpSxevBi77rorAP/xSutRQ+txRTlDLbd5+oyWtxVPs0sb55Qv0JA28u98fOvWrbjllls6/kdl4VEGQZfx5S9/GX19fZgyZQoGBgaGPTqodb2190dKb1POPAMjX/KsvZ+Sd6P5TLcsn7+SjXuZvFvOvVseJ9uV88rdej72ybvqmzZtwne+8x08/vjj+M///M82nIGRhFAGQZdw5ZVXYs8998SkSZOGCRvvgkqBkhMmclLFWopjTbTwcUYu0FnQtOfEebdZyzt3t7n92otbch657rLrLV8C/OCDD+KMM87oiFiGUAZBwyxevBj77LMP9t9//yEvUK5h5IIHlD1KYMcEDx/z4y/WleOMOQ0f6+QCx7vIOb30Drkny8dD+RpPjrSHi2y2WYuf2+LBBx/E1KlTcfbZZ+Pzn//8+J0UQQhlEDTIl770JUyfPn2YN2W9ik0u4eHH5RgiX8Ij0/JwOUHERZJPylhepBRCfizXh3vFMq30IjWPl9sDDJ+Fz/ncdttt6Ovrw+WXXz7GM6ITQhkEDXHVVVdh2rRpw4QCqLcGUXqCXER4upLAafFk117mLUVPprO64dz7lPHlbDovhwuvZe+LL76INWvWmG09VkIog6DD3HTTTejt7cVBBx00bMmOtdRHLrWRfwQm0/GubO5u57QZ7o1p23LiRnbTcxnZDvn+SS0tX4cpu/By0XrOO+cr12fKtunp6cEDDzzQtnMWQhkEHWLx4sXYb7/9sO+++w7zBPn6Q2CHSEgPiguo9mQN99RyOI+nxeFrKmU5sourPYGT4el4Hfhx6TlyrG61fMoo58W/5XY7iHWUQdABrr32Whx99NE44IADzFehcY+Qh2sf7c/B+IeLlgznbwiS6xu1Pxvjdmnla8d5OfL5bisvGScfk39KBgxfp9nX1ze03671leFRBkEbWbx4Md7whjdg+vTpQ15bRnp1fF8bk5NjgNpYn5x4kct8cj7abLX05DLcE815cJtzGaWnd7QhAG1ck4drXqkc6+T1vOeee1o7QTUJoQyCNnH55Zfj7W9/uytEcuwO0IVEdnmtCQ6tey5noGU+0jYurFJouc08vayLZr88bsXX2kl2vflsekppaKz3Xe961yjPlk8IZRC0gaVLl+K4444btuSHT0ho6yL5iyayEOTnofkkR/7OL5vgfwSW8V7mm48DI8cG5RpMjnzqhr8UA4D6Qo3sRfM/OePHNaHlT+xkW+QLNGQ9jjrqKHz3u98tn5hREmOUQdAGTjrpJBx//PFDQiLHDPM373LLD48LjPzDMOmt8Tg8by2OHMuUecv8ZflaOj7WKONotvE20PLPx/MYpxY/77dTJIEQyiBoKy+++KIqgppYaHE0QQN0QZRCpwlySZx5uTltFqtcLrcpY9nQihDL+NIj5vt8ounNb37zeJwqlxDKIGgjM2fOxOrVq4e6pXJWOAuQNbudBcOaXdbi53y5l8ePSxHWxFib9e7r61Nn6+vMzssy+Z+P8bytOmqz8bvvvjve97734aabbmr7eQyhDII2c8EFF2DFihUAbO/I6gbLbqgWLsVPhmtxcpg1NFDqlpfs1brI0iYtrRwysIYUAOCEE07o2H/qxGROEHSAuXPn4utf/zoOOuggAMNfZCGfWCm9WEK+gEJOsmTk8iH5+KBMk8uv8xo27WkdPlnDn7bRFtdzO+QLhXN6/iQOt+/d7343Pve5z43xjLRGeJRB0CHmzJmDX/3qV6pHlfcBe0KGH9Pi1BmLtPKR3e66HqhMz/OXf3amxdHqaZUzadIkHHrooR0XSSCEMgg6yqxZs0Z0RTWhKHVJrW6uFc/qKmuTMnxbirgVxxJP7Z8hZd5evfNn+vTpOP7443HDDTeM9RSMiuh6B0GHOeGEE3DDDTdg2rRpw7rcwKBYyP/ilp/8R118zSPvHvNjvBudj+c0Wtfb65rzbjLvUntdcLkWM6/91NZRyq59SgmTJk3CEUcc0bbXp9UlPMogaIDZs2dj48aNAEZ6UtZSHK0b7XmbVloZx/JSte6v1U3WJoV4nqV9zSM99thjcfvttzcukkAIZRA0xrp168yuplxDqHVrtWPeTHKdOFZ4LsfqJsswazmS3M/wPPr6+jBt2rSuEMhMdL2DoCEWLFiAb37zm5g2bdqILnEWFf5YIEd2t+V7JOWjkzwNnzXPaC+ykI9F8mfF+TPYeZs/E57D5Hstpf28jn/yJ3+CadOm4corrxyfBh5Hih4lEU0lonuIaA0RPUJEF1bhC4noCSJaVX1msjSXENFaInqMiE5pZwWCYCKzYcMGbN261fTGrK6x9OZ4fC1uxkrLP143X3qxOT3/zsc5vCxuL1VDDZMnT8YhhxzSlSIJ1PMo+wFcnFJ6kIj2APAAEd1VHftiSukqHpmI3gTgLABvBvA6AHcT0f9MKQ1/x1QQBJg/fz4A4Atf+AJmzJgxJCh8goO/wYePX8o368h1jPJVZDmetf6Se5DZE+TiKL3MvFYS2LHmUgokXzuZy+evSUspYcqUKfjnf/7ncW7Z8aXoUaaUnkopPVht/xHAGgAHOElOB3BjSunllNIvAawFcOx4GBsEOysXX3wx+vv7h4VJz1LzxqxxTW8c0xqv1NZtyvyk18kfwdTsyMe0vHp7e/Ha176260USaHEyh4imAzgSwPIq6HwieoiIriOiV1dhBwD4NUu2Eb6wBkGAwcmd7LFpogcMX5SurYm01klqXWbZbefxtTcBWd1wy1a+zd8AlMMPPvhgXHPNNe1t1HGi9mQOEe0O4D8AXJRSeoGIvgLgMgCp+v4CgA8D0P68IskAIjoHwDkAMG3atNYtD4KdjI985CMAgHvvvXeYV1f6I7C6jzDmrq/clvnkMN6dl111L03Ov7+/f9iLdfP21KlTsWjRonY25bhTy6Mkol0wKJL/llK6GQBSSk+nlLanlAYAXIsd3euNAKay5K8H8KTMM6X01ZTSjJTSjClTpoylDkGwU7Fx40aza215eXLix1rWk/f5tveopNXd18JlWbLL3dPTgz333HPCiSRQb9abAHwNwJqU0tUsfH8W7UwAD1fbtwI4i4gmE9GBAA4FcP/4mRwEOzfr1q3DCy+80NI4ZEkUNaH1uvF1uvPa44rWbP0uu+yC/fbbD1dffTUmInW63u8EMAvAT4loVRU2D8DZRHQEBrvV6wGcCwAppUeI6LsAfobBGfPzYsY7COqzYMECAMA3vvENTJ8+HcDwWfD89w/Ajtlvbc1kFioeBuzwAvkbevjaTZmn/IsKXq62njPnkfPcddddu3bZT11IPovZBDNmzEgrV65s2owg6Dr+67/+a0iw8j8d5u8sTFnUcjgPyzPpPIyLoMyTx5Ph+blxLa/t27cPW3iew/fcc0985jOfaaz9WoGIHkgpzdCOxZM5QdDFLFu2DMcdd9xQd5Z7lvyPtni4N6nD3wfJxTav1eRP0UgvlYdnLxPY4b3y90dOmjSpqx5BHCvxrHcQdDHz58/H888/P2ysMCPHJjOtPp8tn6xpZQJJG+981atetVOJJBBCGQRdz4YNG4a2tVlqS+DkOkYvjvfiDF4u91SlgOawhQsXtr1NOk10vYOgy7nwwgsBAPfdd9+IbjAwfOJGHgN2vEeSd6+1CSCOfCEGL4dPJvFu/WWXXTYOte1OwqMMggkC/xsJa7mOPFZaXlR3uRBgPxW0s3qRnBDKIJggaOsrAfvRwVbElIdJrHWWmW5YOdNuousdBBOEvL7yhz/84VD3l89CA8P/SoKvo+RvCMqz07krLd8YlLe17nwOy8uF8tuPdnZCKINggnH//ffj6KOPxqRJk4Y9980XefOF38DwF2BkD5ALJF8knsWTvwyYi+6WLVteMQKZCaEMggnG3LlzAWDo7eh84kXzAr2P9hZz/p156aWX8NJLLw15ta80QiiDYIKyYcOGoX9yzN4gsGNMkXeV5Qt0M/ylwNJ75N8f//jHO1m1riMmc4JggjJ//nxs27bNne3W3iikvRzDW6D+29/+tpkKdhHhUQbBBOY973kPvv71r+OQQw4BYP9vt+xK8y52nsDh4USEbdu24dxzz+1ENbqe8CiDYIIzZ84crF+/XvUGrUcQgR1/NZHj83QpJTz99NPNVKgLCaEMgp2AWbNmmY8vcuTich7Glxg9//zzr7iZbY/oegfBTsLxxx+PG264AdOmTRvxFxF53SVfCsTXRAIYel3arFmzmqxGVxIeZRDsRMyePXvERI22n+Ee6LZt27B+/frmjO9iwqMMgp2MX/ziFzj44INH/HkY0fC3kufjW7duxQc+8IEmTe56wqMMgp2MOXPmYMOGDSOWBfGxyBz2m9/8Bo899ljDFnc/4VEGwU7IrFmz8O1vfxuve93rhhae5+8soP39/bH8pybhUQbBTspf//VfD62jzJ+UElauXInTTjsNZ5xxRtMmThhCKINgJ2b58uXYvn07tm7dihUrVuCUU06JZT+jIIQyCHZi5s6diyeeeAI/+clPhl6mEbROjFEGwU5OrIscO+FRBkEQFAihDIIgKBBCGQRBUCCEMgiCoEAIZRAEQYEQyiAIggIhlEEQBAVCKIMgCAqQ/C+NRowgehbAJgDPNW0L4zUIezy6zR6g+2wKe3y6zZ43pJSmaAe6QigBgIhWppRmNG1HJuzx6TZ7gO6zKezx6TZ7PKLrHQRBUCCEMgiCoEA3CeVXmzZAEPb4dJs9QPfZFPb4dJs9Jl0zRhkEQdCtdJNHGQRB0JU0LpREdCoRPUZEa4mokTeLEtF6IvopEa0iopVV2D5EdBcR/bz6fnWbbbiOiJ4hoodZmGoDDfJPVZs9RERHdciehUT0RNVOq4hoJjt2SWXPY0R0ShvsmUpE9xDRGiJ6hIgurMIbaSPHnkbaiIh2JaL7iWh1Zc9nqvADiWh51T7fIaJJVfjkan9tdXz6eNpTsOkbRPRL1kZHVOFt/12PGv4n6Z3+AOgF8AsABwGYBGA1gDc1YMd6AK8RYZ8HMLfangvgyjbbcDyAowA8XLIBwEwASwAQgHcAWN4hexYC+JgS903VuZsM4MDqnPaOsz37Aziq2t4DwONVuY20kWNPI21U1XP3ansXAMuren8XwFlV+L8A+D/V9kcB/Eu1fRaA77ThN2TZ9A0Af6XEb/vverSfpj3KYwGsTSmtSyltBXAjgNMbtilzOoDrq+3rAbT1n5hSSvcB+F1NG04HcEMa5McA9iai/Ttgj8XpAG5MKb2cUvolgLUYPLfjac9TKaUHq+0/AlgD4AA01EaOPRZtbaOqni9Wu7tUnwTg3QBuqsJl++R2uwnASZT/nrH9Nlm0/Xc9WpoWygMA/Jrtb4T/Y2sXCcCdRPQAEZ1The2XUnoKGLwoAOzbgF2WDU222/lVt+g6NhzRUXuqbuKRGPRQGm8jYQ/QUBsRUS8RrQLwDIC7MOi1Pp9S6lfKHLKnOv4HAP9jPO3RbEop5Tb6bNVGXySiydImxd5GaVootTtYE9Pw70wpHQXgNADnEdHxDdjQCk2121cAHAzgCABPAfhCp+0hot0B/AeAi1JKL3hRO2GTYk9jbZRS2p5SOgLA6zHorR7mlNmR9pE2EdFbAFwC4H8BOAbAPgA+2UmbRkPTQrkRwFS2/3oAT3baiJTSk9X3MwBuweCP7Ons9lffz3TaLseGRtotpfR09cMfAHAtdnQdO2IPEe2CQVH6t5TSzVVwY22k2dN0G1U2PA/gXgyO8+1NRPlPBHmZQ/ZUx/dC/aGWsdh0ajVskVJKLwP4Ohpoo1ZpWihXADi0mpmbhMFB5Vs7aQAR7UZEe+RtAO8F8HBlx4eqaB8C8P1O2lVh2XArgNnVLOE7APwhdz/biRgvOhOD7ZTtOauaST0QwKEA7h/nsgnA1wCsSSldzQ410kaWPU21ERFNIaK9q+1XAXgPBsdN7wHwV1U02T653f4KwA9TNaPSZpseZTc2wuCYKW+jjv+ua9H0bBIGZ7oex+B4yqcaKP8gDM5GrgbwSLYBg+M1SwH8vPrep812/DsGu2rbMHhn/RvLBgx2Uf5v1WY/BTCjQ/Z8syrvIQz+qPdn8T9V2fMYgNPaYM9xGOyGPQRgVfWZ2VQbOfY00kYA3grgJ1W5DwNYwH7f92Nw8uh7ACZX4btW+2ur4we14ZxZNv2waqOHAXwLO2bG2/67Hu0nnswJgiAo0HTXOwiCoOsJoQyCICgQQhkEQVAghDIIgqBACGUQBEGBEMogCIICIZRBEAQFQiiDIAgK/H9T06swdEi7/gAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "from matplotlib import pyplot as plt\n",
    "\n",
    "filename = 'pic6.png'\n",
    "origin = cv.imread(filename)\n",
    "gray = cv.cvtColor(origin, cv.COLOR_BGR2GRAY)\n",
    "plt.imshow(gray,'gray')\n",
    "\n",
    "hist = cv.calcHist([gray], [0], None, [256], [0, 256])\n",
    "plt.figure()\n",
    "plt.title(\"Grayscale Histogram\")\n",
    "plt.xlabel(\"Bins\")\n",
    "plt.ylabel(\"# of Pixels\")\n",
    "plt.plot(hist)\n",
    "plt.xlim([0, 256])\n",
    "plt.show()\n",
    "\n",
    "_, otsu = cv.threshold(gray, 0, 255, cv.THRESH_OTSU)\n",
    "plt.imshow(otsu, 'gray')\n",
    "\n",
    "_, flood, _, _ = cv.floodFill(gray, None, (50, 50), 255)\n",
    "plt.imshow(flood, 'gray')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table width=\"100%\">\n",
    "    <tr>\n",
    "        <th width=\"50%\" align=\"center\">灰度图像</th>\n",
    "        <th width=\"50%\" align=\"center\">灰度直方图</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"pic6_gray.png\" /></td>\n",
    "        <td><img src=\"pic6_gray_hist.png\" /></td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <th width=\"50%\" align=\"center\">大津算法</th>\n",
    "        <th width=\"50%\" align=\"center\">区域生长法</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"pic6_gray_otsu.png\" /></td>\n",
    "        <td><img src=\"pic6_gray_flood.png\" /></td>\n",
    "    </tr>\n",
    "</table>\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "4. 使用米粒图像，分割得到各米粒，首先计算各区域(米粒)的面积、长度等信息，进一步计算面积、长度的均值及方差，分析落在3sigma范围内米粒的数量。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "mean of length: 29.651786, variance of length: 49.129460\n",
      "mean of area  : 183.821053, variance of area  : 4367.799557\n",
      "total: 95\n",
      "sigma of length: 7.009241, length in 3 sigma: 94\n",
      "sigma of area: 66.089330, area in 3 sigma: 93\n"
     ]
    }
   ],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "import cmath\n",
    "\n",
    "filename = \"mili.png\"\n",
    "origin = cv.imread(filename)\n",
    "gray = cv.cvtColor(origin, cv.COLOR_BGR2GRAY)\n",
    "cv.imwrite(\"mili_gray.png\", gray);\n",
    "\n",
    "#_, bw = cv.threshold(gray, 0, 255, cv.THRESH_OTSU)\n",
    "gray_blur = cv.GaussianBlur(gray, (3,3), 1)\n",
    "bw = cv.adaptiveThreshold(gray_blur, 255, cv.BORDER_REPLICATE, cv.THRESH_BINARY, 65, 1)\n",
    "cv.imwrite(\"mili_gray_otsu.png\", bw)\n",
    "cv.imshow(\"otsu\", bw)\n",
    "\n",
    "element = cv.getStructuringElement(cv.MORPH_CROSS, (3,3))\n",
    "bw = cv.morphologyEx(bw, cv.MORPH_OPEN, element)\n",
    "\n",
    "seg=bw.copy()\n",
    "\n",
    "_, cnts, _ = cv.findContours(seg, cv.RETR_EXTERNAL, cv.CHAIN_APPROX_SIMPLE)\n",
    "\n",
    "count = 0\n",
    "\n",
    "area_list = []\n",
    "len_list = []\n",
    "\n",
    "for i in range(len(cnts) - 1, -1, -1):\n",
    "    c = cnts[i]\n",
    "    area = cv.contourArea(c)\n",
    "    if (area < 10):\n",
    "        continue\n",
    "        \n",
    "    count = count + 1\n",
    "    #print(\"blob %d: %f\" % (i, area))\n",
    "    \n",
    "    rect = cv.boundingRect(c)\n",
    "    \n",
    "    cv.rectangle(origin, (rect[0], rect[1]), (rect[0]+rect[2], rect[1]+rect[3]), (0,0,0xff), 1)\n",
    "    \n",
    "    strCount = \"%d\" % count\n",
    "    \n",
    "    cv.putText(origin, strCount, (rect[0], rect[1]), cv.FONT_HERSHEY_PLAIN, 0.5, (0, 0xff, 0))\n",
    "    \n",
    "    area_list.append(area)\n",
    "    mili_len = (rect[2]**2+rect[3]**2)**0.5\n",
    "    len_list.append(mili_len)\n",
    "    \n",
    "cv.imshow(\"origin\", origin)\n",
    "cv.imwrite(\"mili_count.png\", origin)\n",
    "\n",
    "len_mean = np.mean(len_list)\n",
    "len_variance = np.var(len_list)\n",
    "area_mean = np.mean(area_list)\n",
    "area_variance = np.var(area_list)\n",
    "\n",
    "print(\"mean of length: %f, variance of length: %f\" % (len_mean, len_variance))\n",
    "print(\"mean of area  : %f, variance of area  : %f\" % (area_mean, area_variance))\n",
    "\n",
    "sigma_len = len_variance ** 0.5\n",
    "sigma_area = area_variance ** 0.5\n",
    "\n",
    "len_3sigma_count = 0\n",
    "area_3sigma_count = 0\n",
    "\n",
    "for l in len_list:\n",
    "    if abs(l - len_mean) < 3*sigma_len:\n",
    "        len_3sigma_count = len_3sigma_count + 1\n",
    "\n",
    "for a in area_list:\n",
    "    if abs(a - area_mean) < 3*sigma_area:\n",
    "        area_3sigma_count = area_3sigma_count + 1\n",
    "\n",
    "print(\"total: %d\" % count)\n",
    "print(\"sigma of length: %f, length in 3 sigma: %d\" % (sigma_len, len_3sigma_count))\n",
    "print(\"sigma of area: %f, area in 3 sigma: %d\" % (sigma_area, area_3sigma_count))\n",
    "\n",
    "key = cv.waitKey()\n",
    "while (key != 32):\n",
    "    key = cv.waitKey()\n",
    "cv.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "大津分割：\n",
    "![米粒图大津分割](mili_gray_otsu.png)\n",
    "统计结果：\n",
    "![米粒图统计结果](mili_count.png)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "扩展作业：\n",
    "\n",
    "5. 使用棋盘格及自选风景图像，分别使用SIFT、FAST及ORB算子检测角点，并比较分析检测结果。  \n",
    "(可选)使用Harris角点检测算子检测棋盘格，并与上述结果比较。 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "\n",
    "def drawOnImage(image, points, radius, color, thickness):\n",
    "    #print(len(points))\n",
    "    for i in points:\n",
    "        #print(i)\n",
    "        cv.circle(image, i, radius, color, thickness)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "\n",
    "filename = \"chessboard.png\"\n",
    "origin = cv.imread(filename)\n",
    "gray = cv.cvtColor(origin, cv.COLOR_BGR2GRAY)\n",
    "cv.imwrite(\"chessboard_gray.png\", gray);\n",
    "\n",
    "cv.imshow(\"origin\", origin)\n",
    "cv.imshow(\"gray\", gray)\n",
    "\n",
    "cornerStrength = cv.cornerHarris(gray, 3, 3, 0.01)\n",
    "_, maxStrength, _, _ = cv.minMaxLoc(cornerStrength)\n",
    "dilated = cv.dilate(cornerStrength, cv.getStructuringElement(cv.MORPH_RECT,(3,3)))\n",
    "localMax = cv.compare(cornerStrength, dilated, cv.CMP_EQ)\n",
    "\n",
    "qualityLevel = 0.5\n",
    "threshold = qualityLevel * maxStrength\n",
    "_, cornerTh = cv.threshold(cornerStrength, threshold, 255, cv.THRESH_BINARY)\n",
    "\n",
    "cornerTh = np.array(cornerTh)\n",
    "cornerMap = cornerTh.astype(np.uint8)\n",
    "result = cornerTh.astype(np.uint8)\n",
    "cv.bitwise_and(cornerMap, localMax, cornerMap)\n",
    "cv.bitwise_and(result, cornerMap, result)#有必要吗？\n",
    "cv.imwrite(\"chessboard_harris.png\", result)\n",
    "cv.imshow(\"harris\", result)\n",
    "\n",
    "points=[]\n",
    "for y in range(0, result.shape[0]):\n",
    "    for x in range(0, result.shape[1]):\n",
    "        if result[y, x] != 0:\n",
    "            points.append((x,y))\n",
    "        \n",
    "drawOnImage(origin, points, 3, (0, 0, 255), 1)\n",
    "cv.imwrite(\"chessboard_harris_mix.png\", origin)\n",
    "cv.imshow(\"result\", origin)\n",
    "\n",
    "key = cv.waitKey()\n",
    "while (key != 32):\n",
    "    key = cv.waitKey() \n",
    "cv.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "棋盘格Harris检测结果：\n",
    "\n",
    "原图：![原图](chessboard.png)\n",
    "Harris角点：![Harris](chessboard_harris.png)\n",
    "合成图：![Mix](chessboard_harris_mix.png)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "\n",
    "filename = \"home.jpg\"\n",
    "origin = cv.imread(filename)\n",
    "gray = cv.cvtColor(origin, cv.COLOR_BGR2GRAY)\n",
    "cv.imwrite(\"home_gray.jpg\", gray);\n",
    "\n",
    "cv.imshow(\"origin\", origin)\n",
    "cv.imshow(\"gray\", gray)\n",
    "\n",
    "cornerStrength = cv.cornerHarris(gray, 3, 3, 0.01)\n",
    "_, maxStrength, _, _ = cv.minMaxLoc(cornerStrength)\n",
    "dilated = cv.dilate(cornerStrength, cv.getStructuringElement(cv.MORPH_RECT,(3,3)))\n",
    "localMax = cv.compare(cornerStrength, dilated, cv.CMP_EQ)\n",
    "\n",
    "qualityLevel = 0.1\n",
    "threshold = qualityLevel * maxStrength\n",
    "_, cornerTh = cv.threshold(cornerStrength, threshold, 255, cv.THRESH_BINARY)\n",
    "\n",
    "cornerTh = np.array(cornerTh)\n",
    "cornerMap = cornerTh.astype(np.uint8)\n",
    "result = cornerTh.astype(np.uint8)\n",
    "cv.bitwise_and(cornerMap, localMax, cornerMap)\n",
    "cv.bitwise_and(result, cornerMap, result)#有必要吗？\n",
    "cv.imwrite(\"home_harris.jpg\", result)\n",
    "cv.imshow(\"harris\", result)\n",
    "\n",
    "points=[]\n",
    "for y in range(0, result.shape[0]):\n",
    "    for x in range(0, result.shape[1]):\n",
    "        if result[y, x] != 0:\n",
    "            points.append((x, y))\n",
    "        \n",
    "drawOnImage(origin, points, 3, (0, 0, 255), 1)\n",
    "cv.imwrite(\"home_harris_mix.jpg\", origin)\n",
    "cv.imshow(\"result\", origin)\n",
    "\n",
    "key = cv.waitKey()\n",
    "while (key != 32):\n",
    "    key = cv.waitKey() \n",
    "cv.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "风景图Harris检测结果：\n",
    "\n",
    "原图：![](home.jpg)\n",
    "Harris角点: ![](home_harris.jpg)\n",
    "合成图：![](home_harris_mix.jpg)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "\n",
    "def detect(feature, image, fname):\n",
    "    keypoints = feature.detect(image)\n",
    "    cv.drawKeypoints(image, keypoints, image, color, cv.DRAW_MATCHES_FLAGS_DRAW_OVER_OUTIMG)\n",
    "    cv.imwrite(fname+\".png\", image)\n",
    "    cv.imshow(fname, image)\n",
    "   \n",
    "def detect_comp(im, imname):\n",
    "    image = im.copy()\n",
    "    feature = cv.ORB_create()\n",
    "    detect(feature, image, imname+\"_ORB\")\n",
    "\n",
    "    image = im.copy()\n",
    "    feature = cv.FastFeatureDetector_create()\n",
    "    detect(feature, image, imname+\"_FAST\")\n",
    "\n",
    "    image = im.copy()\n",
    "    feature = cv.xfeatures2d_SURF.create(1500)\n",
    "    detect(feature, image, imname+\"_SURF\")\n",
    "\n",
    "    image = im.copy()\n",
    "    feature = cv.xfeatures2d_SIFT.create()\n",
    "    detect(feature, image, imname+\"_SIFT\")\n",
    "\n",
    "    image = im.copy()\n",
    "    feature = cv.xfeatures2d_HarrisLaplaceFeatureDetector.create()\n",
    "    detect(feature, image, imname+\"_HLFD\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "\n",
    "filename = \"chessboard.png\"\n",
    "origin = cv.imread(filename)\n",
    "gray = cv.cvtColor(origin, cv.COLOR_BGR2GRAY)\n",
    "cv.imwrite(\"chessboard_gray.png\", gray);\n",
    "\n",
    "cv.imshow(\"origin\", origin)\n",
    "cv.imshow(\"gray\", gray)\n",
    "\n",
    "nfeatures = 128\n",
    "scaleFactor = 2\n",
    "nlevels = 5\n",
    "color = (255, 0, 0)\n",
    "\n",
    "\n",
    "im = origin\n",
    "imname = \"chess_color\"\n",
    "detect_comp(im, imname)\n",
    "\n",
    "im = gray\n",
    "imname = \"chess_gray\"\n",
    "detect_comp(im, imname)\n",
    "\n",
    "key = cv.waitKey()\n",
    "while (key != 32):\n",
    "    key = cv.waitKey() \n",
    "cv.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table width=\"100%\">\n",
    "    <tr>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格彩色图像</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格彩色ORB检测</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格彩色FAST检测</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"chessboard.png\" /></td>\n",
    "        <td><img src=\"chess_color_ORB.png\" /></td>\n",
    "        <td><img src=\"chess_color_FAST.png\" /></td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格彩色SURF检测</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格彩色SIFT检测</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格彩色HLFD检测</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"chess_color_SURF.png\" /></td>\n",
    "        <td><img src=\"chess_color_SIFT.png\" /></td>\n",
    "        <td><img src=\"chess_color_HLFD.png\" /></td>\n",
    "    </tr>\n",
    "</table>\n",
    "<table width=\"100%\">\n",
    "    <tr>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格灰度图像</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格灰度ORB检测</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格灰度FAST检测</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"chessboard_gray.png\" /></td>\n",
    "        <td><img src=\"chess_gray_ORB.png\" /></td>\n",
    "        <td><img src=\"chess_gray_FAST.png\" /></td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格灰度SURF检测</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格灰度SIFT检测</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格灰度HLFD检测</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"chess_gray_SURF.png\" /></td>\n",
    "        <td><img src=\"chess_gray_SIFT.png\" /></td>\n",
    "        <td><img src=\"chess_gray_HLFD.png\" /></td>\n",
    "    </tr>\n",
    "</table>\n",
    "\n",
    "HLFD应用于棋盘格灰度图像未检测到任何特征点"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "\n",
    "filename = \"home.jpg\"\n",
    "origin = cv.imread(filename)\n",
    "gray = cv.cvtColor(origin, cv.COLOR_BGR2GRAY)\n",
    "cv.imwrite(\"home_gray.png\", gray);\n",
    "\n",
    "cv.imshow(\"origin\", origin)\n",
    "cv.imshow(\"gray\", gray)\n",
    "\n",
    "nfeatures = 128\n",
    "scaleFactor = 2\n",
    "nlevels = 5\n",
    "color = (255, 0, 0)\n",
    "\n",
    "im = origin\n",
    "imname = \"home_color\"\n",
    "detect_comp(im, imname)\n",
    "\n",
    "im = gray\n",
    "imname = \"home_gray\"\n",
    "detect_comp(im, imname)\n",
    "\n",
    "key = cv.waitKey()\n",
    "while (key != 32):\n",
    "    key = cv.waitKey() \n",
    "cv.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<table width=\"100%\">\n",
    "    <tr>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格彩色图像</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格彩色ORB检测</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格彩色FAST检测</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"home.jpg\" /></td>\n",
    "        <td><img src=\"home_color_ORB.png\" /></td>\n",
    "        <td><img src=\"home_color_FAST.png\" /></td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格彩色SURF检测</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格彩色SIFT检测</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格彩色HLFD检测</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"home_color_SURF.png\" /></td>\n",
    "        <td><img src=\"home_color_SIFT.png\" /></td>\n",
    "        <td><img src=\"home_color_HLFD.png\" /></td>\n",
    "    </tr>\n",
    "</table>\n",
    "<table width=\"100%\">\n",
    "    <tr>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格灰度图像</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格灰度ORB检测</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格灰度FAST检测</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"home_gray.jpg\" /></td>\n",
    "        <td><img src=\"home_gray_ORB.png\" /></td>\n",
    "        <td><img src=\"home_gray_FAST.png\" /></td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格灰度SURF检测</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格灰度SIFT检测</th>\n",
    "        <th width=\"33%\" align=\"center\">棋盘格灰度HLFD检测</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"home_gray_SURF.png\" /></td>\n",
    "        <td><img src=\"home_gray_SIFT.png\" /></td>\n",
    "        <td><img src=\"home_gray_HLFD.png\" /></td>\n",
    "    </tr>\n",
    "</table>\n",
    "\n",
    "从图片看，FAST算法得到的特征点最多，其次是SURF、SIFT、HLFD，最好的是ORB，保留了前景主要特征点又排除了天空中的飞鸟和云朵\n",
    "\n",
    "对棋盘格灰度图像不起作用的HLFD算法在这个图片的灰度图像还是有不错的效果的，对彩色图像效果不理想"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "\n",
    "filename = \"home.jpg\"\n",
    "origin = cv.imread(filename)\n",
    "\n",
    "def detect(hessianThreshold, im, fname):\n",
    "    image = im.copy()\n",
    "    feature = cv.xfeatures2d_SURF.create(hessianThreshold)\n",
    "    keypoints = feature.detect(image)\n",
    "    cv.drawKeypoints(image, keypoints, image, (255, 0, 0), cv.DRAW_MATCHES_FLAGS_DRAW_OVER_OUTIMG)\n",
    "    cv.imwrite(fname+\".jpg\", image)\n",
    "    cv.imshow(fname, image)\n",
    "    \n",
    "detect(500, origin, \"home_SURF_500\")\n",
    "detect(1000, origin, \"home_SURF_1000\")\n",
    "detect(1500, origin, \"home_SURF_1500\")\n",
    "detect(2000, origin, \"home_SURF_2000\")\n",
    "detect(2500, origin, \"home_SURF_2500\")\n",
    "detect(3000, origin, \"home_SURF_3000\")\n",
    "    \n",
    "key = cv.waitKey()\n",
    "while (key != 32):\n",
    "    key = cv.waitKey() \n",
    "cv.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "SURF算法中不同的hessianThreshold值对比：\n",
    "<table width=\"100%\">\n",
    "    <tr>\n",
    "        <th width=\"50%\" align=\"center\">原图</th>\n",
    "        <th width=\"50%\" align=\"center\">hessianThreshold=100</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"home.jpg\" /></td>\n",
    "        <td><img src=\"home_color_SURF.png\" /></td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <th width=\"50%\" align=\"center\">hessianThreshold=500</th>\n",
    "        <th width=\"50%\" align=\"center\">hessianThreshold=1000</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"home_SURF_500.jpg\" /></td>\n",
    "        <td><img src=\"home_SURF_1000.jpg\" /></td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <th width=\"50%\" align=\"center\">hessianThreshold=1500</th>\n",
    "        <th width=\"50%\" align=\"center\">hessianThreshold=2000</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"home_SURF_1500.jpg\" /></td>\n",
    "        <td><img src=\"home_SURF_2000.jpg\" /></td>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <th width=\"50%\" align=\"center\">hessianThreshold=2500</th>\n",
    "        <th width=\"50%\" align=\"center\">hessianThreshold=3000</th>\n",
    "    </tr>\n",
    "    <tr>\n",
    "        <td><img src=\"home_SURF_2500.jpg\" /></td>\n",
    "        <td><img src=\"home_SURF_3000.jpg\" /></td>\n",
    "    </tr>\n",
    "</table>\n",
    "\n",
    "hessianThreshold值越大，保留的特征点越少，"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "keyPoint1.size =  1250\n",
      "keyPoint2.size =  1226\n",
      "keyPoint1.size =  1104\n",
      "keyPoint2.size =  1086\n",
      "keyPoint1.size =  1046\n",
      "keyPoint2.size =  1000\n",
      "keyPoint1.size =  50\n",
      "keyPoint2.size =  50\n",
      "keyPoint1.size =  50\n",
      "keyPoint2.size =  50\n"
     ]
    }
   ],
   "source": [
    "import cv2 as cv\n",
    "import numpy as np\n",
    "\n",
    "fn1 = \"left03.jpg\"\n",
    "fn2 = \"left14.jpg\"\n",
    "im1 = cv.imread(fn1)\n",
    "im2 = cv.imread(fn2)\n",
    "\n",
    "def match(feature, im1, im2, outname):\n",
    "    keyPoint1 = feature.detect(im1)\n",
    "    keyPoint2 = feature.detect(im2)\n",
    "    _, descriptors1 = feature.compute(im1, keyPoint1)\n",
    "    _, descriptors2 = feature.compute(im2, keyPoint2)\n",
    "    matcher = cv.DescriptorMatcher_create(\"BruteForce\")\n",
    "    matches = matcher.match(descriptors1, descriptors2)\n",
    "    img_matches = np.empty(im2.shape)\n",
    "    img_matches = cv.drawMatches(im1, keyPoint1, im2, keyPoint2, matches, img_matches)\n",
    "    cv.imwrite(outname+\".png\", img_matches)\n",
    "    cv.imshow(outname, img_matches)\n",
    "    print(\"keyPoint1.size = \", len(keyPoint1))\n",
    "    print(\"keyPoint2.size = \", len(keyPoint2))\n",
    "    \n",
    "feature = cv.xfeatures2d_SURF.create(1000)\n",
    "match(feature, im1, im2, \"chess_match_SURF_1000\")\n",
    "feature = cv.xfeatures2d_SURF.create(2000)\n",
    "match(feature, im1, im2, \"chess_match_SURF_2000\")\n",
    "feature = cv.xfeatures2d_SURF.create(3000)\n",
    "match(feature, im1, im2, \"chess_match_SURF_3000\")\n",
    "feature = cv.ORB_create(50)\n",
    "match(feature, im1, im2, \"chess_match_ORB_50\")\n",
    "\n",
    "fn1 = \"ela_00.jpg\"\n",
    "fn2 = \"ela_90.jpg\"\n",
    "im1 = cv.imread(fn1)\n",
    "im2 = cv.imread(fn2)\n",
    "match(feature, im1, im2, \"ela_match_ORB_50\")\n",
    "\n",
    "key = cv.waitKey()\n",
    "while (key != 32):\n",
    "    key = cv.waitKey() \n",
    "cv.destroyAllWindows()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![ela_ORB match](ela_match_ORB_50.png)\n",
    "![chess ORB match](chess_match_ORB_50.png)\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
